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ABSTRACT 

We present the first results of hydrodynamical simulations that follow the formation of galax- 
ies to the present day in nearly spherical regions of radius 2Qh~^ Mpc drawn from the 
Millennium Simulation (Springel et al.). The regions have mean overdensities that deviate by 
(—2, —1, 0, +1,+2)(T from the cosmic mean, where a is the rms mass fluctuation on a scale 
of ^ 20 /i^^ Mpc at z = 1.5. The simulations have mass resolution of up to ^ 10*' hr"^ Mq, 
cover the entire range of large-scale cosmological environments, including rare objects such as 
massive clusters and sparse voids, and allow extrapolation of statistics to the (500 h^^ Mpc)^ 
Millennium Simulation volume as a whole. They include gas cooling, photoheating from an 
imposed ionising background, supemova feedback and galactic winds, but no AGN. In this 
paper we focus on the star formation properties of the model. We find that the specific star for- 
mation rate density at z < 10 varies systematically from region to region by up to an order of 
magnitude, but the global value, averaged over all volumes, closely reproduces observational 
data. Massive, compact galaxies, similar to those observed in the GOODS fields (Wiklind et 
al.), form in the overdense regions as early as z = 6, but do not appear in the underdense 
regions until z ^ 3. These environmental variations are not caused by a dependence of the 
star formation properties on environment, but rather by a strong variation of the halo mass 
function from one environment to another, with more massive haloes forming preferentially 
in the denser regions. At all epochs, stars form most efficiently in haloes of circular velocity 
Vc ~ 250 km . However, the star-formation history exhibits a form of "downsizing" (even 
in the absence of AGN feedback): the stars comprising massive galaxies at z = have mostly 
formed by z = 1 — 2, whilst those comprising smaller galaxies typically form at later times. 
However, additional feedback is required to limit star formation in massive galaxies at late 
times. 

Key words: galaxies: abundances - galaxies: clusters: general - galaxies: formation - galax- 
ies: intergalactic medium - methods: N-hoAy simulations 



1 INTRODUCTION 
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Numerical simulations have emerged, over the past two decades or 
so, as a useful technique for modelling the formation and evolu- 
tion of cosmic structures. In particular, they have yielded accurate 
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predictions for tlie clustering and evolution of dark matter, a rela- 
tively simple problem in which the dynamics are determined purely 
by gravitational forces (e.g. Davis et al. 1985). Modelling the evo- 
lution of small-scale structures involving baryons is considerably 
more challenging because of the complexity of the physics involved 
and the greater uncertainties in the modelling techniques (e.g. Katz 
& Gunn 1991). In the regime where baryonic processes are non- 
negligible, the ability of numerical simulations to yield unique and 
robust predictions is therefore reduced. This is particularly true for 
systems in which radiative cooling is efficient. 

The inherent difficulties and computational expense of mod- 
elling baryons in this regime has stimulated the development of 
semi-analytic techniques (White & Frenk 1991; Kauffmann et al. 
1993; Cole et al. 1994, 2000; Somerville & Primack 1999; Baugh 
2006). These employ simplified prescriptions for the evolution 
of baryons in dark matter haloes. The formation histories of the 
haloes are followed using merger trees constructed anal3^ically 
with Monte Carlo techniques, or extracted directly from A^-body 
simulations. Notable examples of the latter procedure are the stud- 
ies of Bower et al. (2006), Croton et al. (2006), De Lucia et al. 
(2006) and Font et al. (2008), based on merger histories de- 
rived from the Millennium Simulation (Springel et al. 2005). The 
scale and resolution of this dark matter simulation permits semi- 
analytical calculations of the evolution of galaxies more massive 
than the Small Magellanic Cloud within a comoving volimie sim- 
ilar to that probed by the 2dFGRS (CoUess et al. 2001) and SDSS 
(York et al. 2000) at their median redshifts. 

Simulations employing semi-analytic models have yielded 
critical insights into the effects of baryonic processes on the observ- 
able Universe. Nevertheless, the approximations inherent in this 
technique limit the kind of processes that can be followed reliably 
in a relatively simple manner. For example, while semi-analytic 
techniques can be extended to study the evolution of the inter- 
galactic medium (IGM; e.g. Benson et al. 2002), they are not well 
suited for studying the interaction of gas with galactic ejecta, such 
as winds, and the associated dynamical evolution. Processes of this 
kind can only be studied in detail with simulations that explicitly 
follow the full hydrodynamic evolution of the baryons (Wiersma 
et al. 2009b). 

The introduction and subsequent development of hydrody- 
namic numerical methods within the framework of cosmological 
simulations has led to important advances in our understanding 
of structure formation. For example, the relationship between the 
Lyman-a forest and the cosmic large-scale structure was first con- 
vincingly demonstrated using such simulations (e.g. Cen et al. 
1994; Zhang et al. 1995; Hernquist et al. 1996; Miralda-Escude 
et al. 1996; Zhang et al. 1997; Theuns et al. 1998; Dave et al. 
1999). At the opposite end of the density scale, hydrodynamic sim- 
ulations of the formation of individual galaxies have confirmed the 
generic requirement in hierarchical models (first identified using 
semi-analytic techniques) for energy feedback mechanisms to pre- 
vent the over-cooling of gas in small haloes at early times (White 
& Rees 1978; Cole 1991; White & Frenk 1991), and the associ- 
ated transfer of angular momentum from baryons to dark matter 
(e.g. Navarro & Benz 1991; Weil et al. 1998; Thacker & Couchman 
2001; Abadi et al. 2003a,b; Sommer-Larsen et al. 2003; Govemato 
et al. 2004; Robertson et al. 2004; Okamoto et al. 2005; Govemato 
et al. 2007; Scannapieco et al. 2008; Zavala et al. 2008). 

The traditional approach of hydrodynamical simulations has 
involved a trade-off between resolution and computational volume. 
Such trade-offs, however, are not possible if an unbiased statis- 
tical sample of well resolved simulated galaxies is desired. This 



presents a considerable computational challenge, and consequently 
previous simulations have been limited to tracing volumes that are 
relatively small in comparison with those used in large-scale struc- 
ture calculations and those probed by observational surveys. The 
gulf in mass and spatial scale between individual galaxies and 
their large-scale environment compounds the more fundamental 
dynamic range problem faced by simulations of galaxy formation, 
that stems from the influence of individual stars (via stellar winds 
and supernovae) on gas. 

Recent attempts to produce such samples include the 
smoothed particle hydrodynamics (SPH) simulations of Pearce 
et al. (1999, 2001) using the Hydra code (Couchman et al. 1995); 
those of Dave et al. (2006); Oppenheimer & Dave (2006); Dave & 
Oppenheimer (2007); Oppenheimer & Dave (2008a); Croft et al. 
(2008); Oppenheimer & Dave (2008b), using variants of the Gad- 
GEt2 code (Springel 2005); and those of Brooks et al. (2007) using 
the Gasoline code (Wadsley et al. 2004). The adaptive mesh re- 
finement (AMR) technique for hydrodynamics has also been em- 
ployed for simulations of this type, for example by Nagamine 
et al. (2005, see references therein for code details), and by Ocvirk 
et al. (2008), using the Ramses code (Teyssier 2002). To high- 
light an example. Croft et al. (2008) simulated (with 486^ parti- 
cles of gas and dark matter) a periodic volume in which haloes of 
Milky Way type galaxies were resolved with ~ 10^ particles. Al- 
though ~ 1000 galaxies formed in this simulation, its relatively 
small size (a cube of side L = 33.75 Mpc) yielded few mas- 
sive galaxies and sampled a relatively narrow range of cosmologi- 
cal environments. A further technical limitation of relatively small 
volumes is that fluctuations on scales comparable to the size of the 
box become non-negUgible at a relatively high redshift (Bagla & 
Ray 2005; Sirko 2005), after which the simulation can no longer be 
considered a faithful representation of the underlying model. For 
this reason, Croft et al. (2008) halted their simulation at z = 1. 

The need for detailed predictions of the properties of the 
low-redshift IGM and the evolution of galaxies over the interval 
< z < 1 is highlighted by a number of observational develop- 
ments. The fate of the baryons seen in the Lyman-a forest at high 
redshift is still largely uncharted but measurements of the Lyman- 
a forest are beginning to probe the diffuse baryons associated with 
the filaments of the cosmic web (e.g. Pichon et al. 2001; Caucci 
et al. 2008), where a significant fraction may be in the form of a 
"warm-hot intergalactic medium" (WHIM) phase (e.g. Cen & Os- 
triker 1999; Dave et al. 1999). The recent installation of the Cosmic 
Origins Spectrograph (Green et al. 2003) aboard the Hubble Space 
Telescope (HST), may yield significant advances in the understand- 
ing of the low-redshift IGM. Alongside studies of the evolution 
of cosmic gas, realistic modelling of the low-redshift evolution of 
galaxies is essential to coimect the snapshots of cosmic evolution 
provided by galaxy surveys at different epochs, for example DEEP2 
(Davis et al. 2003), VVDS (Le Fevre et al. 2005) and COMBO-17 
(Wolf et al. 2003) at a ~ 1 with the 2dFGRS and SDSS at z - 0. 

Simulating large cosmological volumes (L > lQQh~^ Mpc) 
at high resolution (mgas < 10^ M©) to a = is infeasible 
with current computational resources. To circumvent this funda- 
mental limitation, we have devised a new approach consisting of 
simulating at high resolution regions extracted from the dark matter 
Milleimium Simulation whose overdensities at a = 1.5 represent 
(—2, —1, 0, -l-l,-|-2)iT deviations from the cosmic mean, where a 
is the rms mass fluctuation on a scale of ~ 2Qh^^ Mpc. These 
simulations, the Galaxies-lntergalactic Medium Interaction Calcu- 
lation (GiMiC), foUow the evolution of five roughly spherical re- 
gions of radius ~ 20 Mpc from "zoomed" initial conditions 
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(Frenk et al. 1996; Power et al. 2003; Navarro et al. 2004) using 
full gas dynamics, wiiile tlie rest of tlie Millennium volume is res- 
imulated at lower resolution following dark matter only. We thus 
follow a very wide range of cosmological environments whilst only 
simulating 0.13 per cent of the Millennium volume at high resolu- 
tion. The extent of the dynamic range allowed by our procedure is 
illustrated in Fig. 1 which, starting from the full Millennium Sim- 
ulation volume (L = 500 Mpc), zooms in, first by a factor of 
10 to show a fuU Gimic sphere in the central panel, and then by a 
further factor of 1000 to show a disc galaxy within this region. 

In contrast to simulations of small periodic volumes, our 
adopted strategy allows us to follow the cosmic evolution to z = 0, 
since fluctuations on the scale of the Millennium Simulation are 
still well described by linear theory. Judicious selection of the 
five Gimic regions across enviroimientally diverse regions yields 
a sample spanning the range of structures found within the Millen- 
nium Simulation. In our highest resolution realisations, the Jeans 
scale in the IGM after the epoch of reionisation is marginally re- 
solved. A further advantage of resimulating regions of the Millen- 
nium Simulation is that they complement the existing semi-analytic 
calculations implemented on it (Bower et al. 2006; Croton et al. 
2006; De Lucia et al. 2006; Font et al. 2008). 

The simulations were performed using GADGETS, a sig- 
nificantly upgraded variant of the Gadget2 code described by 
Springel (2005), and including new treatments of baryonic pro- 
cesses, such as radiative cooling, heating by the metagalactic ul- 
traviolet background radiation, star formation, stellar feedback and 
chemodjmamics. The simulations are computationally expensive, 
and so we have limited ourselves to running them with a single im- 
plementation of the code. Analyses of these simulations will there- 
fore, in general, focus on the environmental effects that these sim- 
ulations are able to explore; we do not examine the role of, for 
example, changes to the assumed initial mass function (IMF) or 
the implementation of feedback. Such tests will be explored in the 
closely related Overwhelmingly Large Simulations (OWLS; Schaye 
et al., in prep) project. In this, the first paper in a series, we focus 
on one central aspect of these simulations, the difference in the star 
formation rate density amongst the five different enviroimients, and 
on how these combine to produce a "cosmic" star formation rate 
density (e.g. Lilly et al. 1996; Madau et al. 1996). 

The paper is laid out as follows. In Section 2 we describe the 
simulations, concentrating on new aspects of our code, and provide 
a brief overview of the generation of the initial conditions. Our suite 
of simulations includes runs with varying resolution to enable us to 
check for numerical convergence. In Section 3, we discuss the halo 
mass functions of the five GiMIC regions, which are very different 
from each other. In Section 4, we calculate the evolution of the star 
formation rate of each region and present an estimate for the entire 
Milleimium Simulation based on a weighted average of these mea- 
surements. The star formation properties of haloes and their galax- 
ies are discussed in Section 5. In Section 6, we carry out a limited 
comparison with observational data. We summarise and conclude 
in Section 7. 

The simulations adopt the same cosmological parameters as 
the Millennium Simulation: Om = 0.25, = 0.75, Vlh = 0.045, 
= 1, = 0.9, Ho = 100 h kms"^ Mpc"\ h = 0.73. 
This work is part of the programme of the Virgo consortium for 
cosmological simulations. 



2 METHODOLOGY 

This section describes the three central aspects of our methodology: 
the code used to perform the simulations, the setup of the initial 
conditions, and the identification of dark matter haloes and galax- 
ies. Technical details of the generation of the initial conditions are 
deferred to an Appendix. 

2.1 Simulation code 

For this study, we use GADGET3, an updated version of the Gad- 
GET2 code (Springel 2005), to which several physics modules have 
been added. The domain decomposition differs from that in Gad- 
GEt2 in that it improves load-balancing particularly for simula- 
tions (such as those described here) with strongly clustered particle 
distributions run on parallel supercomputers with a large number 
of cores. The hydrodynamics implementation, taken from Gad- 
Get2, is the entropy conserving formulation of SPH (Gingold & 
Monaghan 1977; Lucy 1977), as discussed in Springel & Hernquist 
(2002). SPH represents a fluid by a set of particles that carry along 
a number of internal properties (for example, mass, entropy), while 
other properties are computed by interpolation over neighbouring 
particles (for example, the density or pressure gradient). 

We have used new physics modules for star formation (Schaye 
& Dalla Vecchia 2008), stellar feedback (Dalla Vecchia & Schaye 
2008), radiative cooling (Wiersma et al. 2009a), and chemodynam- 
ics (Wiersma et al. 2009b). Here we only provide a brief overview 
of the main features of these modules. 

• Gas cooling and photoionisation. Radiative cooling was im- 
plemented as described in Wiersma et al. (2009a)^. In brief, we 
assume the IGM is ionised and heated by a pervading uniform, but 
redshift dependent, ionising background from galaxies and quasars, 
as computed by Haardt & Madau (2001). We assume hydrogen 
reionises at redshift z = and helium II at z = 3.5 (Schaye 
et al. 2000; Theuns et al. 2002). The cooling rate is computed as a 
function of redshift, gas density, temperature and composition on 
an element-by-element basis using interpolation tables computed 
with CLOUDY (as described in Ferland et al. 1998), assuming the 
gas to be optically thin and in ionisation equilibrium. During the 
reionisation of helium II, we inject 2 eV per proton, smoothed with 
a Gaussian function, centred at a = 3.5 and with width (t{z) = 0.5, 
to mimic non-equilibrium and radiative-transfer effects (e.g. Abel 
& Haehnelt 1999). As demonstrated in Wiersma et al. (2009b), with 
this extra heat input the predicted temperature at the mean density is 
in good agreement with the measurements of Schaye et al. (2000). 

• Quiescent star formation and feedback. Cosmological simu- 
lations cannot, at present, resolve the multiphase structure of star- 
forming gas that arises from the combination of gas cooling and 
heating due to massive stars and supernovae (SNe). These pro- 
cesses make star-forming gas resistant to compression, since the en- 
ergy output from SNe rapidly heats the interstellar medium (ISM), 
increasing its pressure. We mimic this by imposing an effective 
equation of state P — Kp^^°s on gas with density exceeding 
riH = 0.1 cm^'\ the critical density for the onset of the thermo- 
gravitational instability (Schaye 2004). Schaye & Dalla Vecchia 
(2008) demonstrate that the exponent 7eos = 4/3 makes both 
the Jeans mass and the ratio of the SPH kernel and the Jeans length 
independent of the gas density, thereby preventing spurious frag- 
mentation due to a lack of numerical resolution. We adopt this 

^ We used their equation (3), rather than (4), and cloudy version 05.07. 
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value here. The equation of state is normahsed such that P/k = 
10^ cm^ ' K for atomic gas at the density threshold. 

Observations of galaxies reveal a tight relationship between gas 
column density and star formation rate, the Kennicutt-Schmidt law 
(e.g. Kennicutt 1998). Schaye & Dalla Vecchia (2008) show that 
since gas column density in a self-gravitating system is directly 
related to pressure, it is possible to implement a local Kennicutt- 
Schmidt law as a pressure law, independently of the chosen value of 
7EOS and even for very low numerical resolution. The star forma- 
tion algorithm therefore converts gas particles with densities above 
the threshold nn = 0.1 cm~^ into stars at a rate that depends on 
the gas pressure, 

m. = A{1 Mq pc-')-'Vrig [-^hP) , (1) 

where mg is the mass of the gas particle for which we are com- 
puting rh*, 7 = 5/3 for a monatomic gas (and should not be con- 
fused with 7eos), /g is the mass fraction in gas (assumed here to be 
unity), and P is the gas pressure. The normalisation. A, and slope, 
n, follow from the observed Kennicutt (1998) law, 

/ y Xl.4 

= 1.5 X 10-^ Mq yr-1 kpc'^ , . " , , (2) 

where we have divided Kennicutt's normalisation by a factor 1.65 
to account for the fact that we assvmie the IMF takes the form pro- 
posed by Chabrier (2003) rather than the commonly adopted form 
due to Salpeter (1955). The conversion of gas particles into stars, 
based upon their associated star formation rate, is stochastic, and 
since the spawning of multiple star particles from a gas particle 
can lead to an artificial reduction in the efficiency of feedback pre- 
scriptions (for a discussion, see Schaye & Dalla Vecchia 2008), we 
convert entire gas particles into stars. Conveniently, this practice 
also conserves the particle number throughout our simulations. 

• Kinetic feedback. Feedback from star formation in the disc 
is represented in part by the imposed relation P = Kp^^°^. 
However, observed starburst galaxies exhibit signs of galaxy-wide 
winds (Heckman et al. 1990; Martin 1999; Pettini et al. 2002; Adel- 
berger et al. 2003; Shapley et al. 2003; Wilman et al. 2005; S win- 
bank et al. 2007) that are thought to be responsible for enriching the 
IGM with metals. We model the generation of winds as follows. Af- 
ter a delay of 3 x 10^ yr, corresponding to the maximum lifetime 
of stars that undergo core collapse SNe, newly formed star parti- 
cles impart a randomly directed 600 km kick to, on average, 

/rh* = 4 of its neighbours. This mass loading value 
was chosen to scale the global star formation rate density, p*{z), 
such that it reasonably matches observational data. 

Assuming that each star with initial mass in the range 6 — 
100 M© injects a kinetic energy of 10®^ erg, these parameter val- 
ues imply that the total wind energy accounts for 80 per cent of the 
available kinetic energy for our chosen IMF. Note that contrary to 
the widely used kinetic feedback recipe of Springel & Hemquist 
(2003a), wind particles are injected locally and are not temporar- 
ily decoupled hydrodynamically. As discussed in Dalla Vecchia & 
Schaye (2008), these changes have a large effect on the structure of 
galaxy discs and outflows. 

• Chemodynamics. Each star particle represents a single stel- 
lar population, with metal abundances inherited from its parent gas 
particle. Given our assumed initial mass function (Chabrier 2003, 
stellar mass range 0.1 — 100 M(:j) and stellar evolution tracks de- 
pendent on metal abundance (Portinari et al. 1998; Mango 2001; 
Thielemann et al. 2003), we follow the delayed release of 11 el- 
ements (hydrogen, helium, carbon, nitrogen, oxygen, neon, mag- 



nesium, silicon, sulphur, calcium, and iron) by type la and type 11 
SNe, and AGB stars. Star particles distribute the synthesised ele- 
ments and lost mass during each timestep to neighbouring gas par- 
ticles using the SPH interpolation scheme. Throughout this study, 
we use the solar abundances of CLOUDY; solar metallicity is there- 
fore taken to be Zq = 0.0127. 



2.2 Initial conditions and run details 

For brevity, we present here only an outline of the generation of the 
initial conditions; a detailed technical description is deferred to the 
Appendix. There, we also discuss how we address the existence of 
an artificial boundary between the high-resolution region and the 
external, low-resolution pressureless region, and how we combine 
results from the five individual GiMic regions to construct esti- 
mates of properties for the entire Millennium Simulation volvune. 

We follow the evolution of a wide range of environments by 
simulating five spherical regions whose overdensities at z = 1.5 
are (—2, —1, 0, -|-1, -|-2) times the root-mean-square deviation, a, 
from the mean on some spatial scale. Ideally the simulations would 
resolve the Jeans mass at our imposed star formation threshold den- 
sity, but this remains beyond the scope of our available computa- 
tional resources. Therefore, for our highest resolution simulation, 
we aimed to resolve the Jeans mass in the post-reionisation IGM, 
which was shown by Theuns et al. (1998) to be a prerequisite for 
attaining converged properties of the Lya forest. This requires a 
gas particle mass of - 10'' ir^ Mq, limiting the size of the re- 
gions to a radius of 18 Mpc. To ensure we also followed a rich 
cluster, we imposed a precondition that the -\-2cr region be centred 
on a massive dark matter halo at a = and increased its radius to 
25/i"^Mpc. 

In order to preserve the correct large-scale forces, the volume 
external to these five regions was simulated at lower resolution and 
without baryonic physics. This methodology is shown schemati- 
cally in Fig. 1 . The multiresolution particle distribution generated 
for each simulation is explained in the Appendix; the particle dis- 
placement fields were calculated using the techniques described in 
Power et al. (2003) and Navarro et al. (2004). 

We created two realisations of the 18/i~^Mpc spheres: 
high-resolution (marginally resolving the IGM Jeans mass) and 
intermediate-resolution (with a factor of 8 fewer particles); we re- 
serve the term "low-resolution" for the original Millennium Sim- 
ulation realisation. The +2a region was generated at intermediate 
resolution only because of its far greater computational demands. 
We have hence carried out 9 simulations in total, with the charac- 
teristics Usted in Table 1. All intermediate resolution realisations 
have been run to = 0. The least cpu-demanding region (— 2cr) 
was also run to « = at high resolution, whilst the three regions, 
(—1,0, +1)(T, were run to z = 2 at high resolution. Simulations of 
the same region at different resolution are needed in order to check 
numerical convergence. In what follows, we perform such tests at 
2 ^ 2 for the (—1, 0, -\-l)a and all the way to z = for the — 2(t 
region. 

In all simulations, the gravitational forces of the baryonic and 
high resolution dark matter particles were softened over the same 
length scale. The softening length is initially fixed in comoving 
space, but becomes fixed in physical space at a predefined redshift, 
i.e. ecom(<i) = rnin(ecom, e™i^s/a)- The softenings were chosen 
such that at « = 3, tiiey are fixed at e^^ys = (1-0, 0.5) kpc for 
the intermediate- and high-resolution runs, respectively. The gas 
particle masses are 1.16 x lO'' M© and 1.45 x 10*^ M© 
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Figure 1. Schematic representation of GiMIC, illustrating the large dynamic range allowed by the "zoomed" initial conditions technique. The left panel shows 
the density distribution within the L = 500 Mpc periodic volume of the Millennium Simulation at z = 1.5. Shown embedded within this volume, to 
scale, is the Oct Gimic region zoomed down to L ~ 50 Mpc and showing the gas density. The right-hand panel, zoomed yet further to L ~ 50 h^^ kpc, 
depicts the gas density within a disc galaxy of mass similar to that of the Milky Way. 
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Table 1. Parameters for the five GiMIC regions. Columns 2-5 give the location (in Millennium Simulation coordinates) and the nominal comoving radius of the 
regions at 2 = 1.5. The following two columns show the number of gas (or, equivalently, dark matter) particles within the zoomed region of the simulation, 
whilst the final column gives the redshift at which each simulation was terminated. Note that the +2a region was not run at high resolution. 

at intermediate and high resolution respectively. These are much 
smaller than the limit required to avoid artificial two-body heat- 
ing effects (Steinmetz & White 1997) that dominate over radiative 
cooling. The dark particle masses are a factor (Sim — ^lh)/^b ~ 
4.56 larger that the gas particle masses. 



2.3 Halo and galaxy identification 

We identify haloes by applying the friends-of-friends (FoF) algo- 
rithm to dark matter particles, using the standard value of the link- 
ing length in units of the mean interparticle separation (b — 0.2; 
Davis et al. 1985). In order to identify the baryonic content of dark 
matter haloes, our group finding algorithm locates the nearest dark 
matter particle to each baryonic (i.e. gas or star) particle. If this dark 
matter particle has been grouped by FoF, the corresponding bary- 
onic particle is also associated with the FoF group. The assignment 
of baryons to dark matter haloes is therefore unambiguous. This 
method occasionally introduces artefacts, mostly for groups with 
few particles. For example, we find that the baryon fraction of low- 
mass haloes varies widely, from haloes with almost no baryons, to 
haloes with more than the cosmic mean value, because our scheme 
allows the baryon distribution around a halo to be more extended 
than the dark matter. Since these problems are mostly restricted to 
small haloes close to the resolution limit, we have not attempted to 
overcome them. 

The FoF algorithm identifies isodensity contours of 5 ~ 



3/(27r&^) ~ 60 (Lacey & Cole 1994), but noise in the particle 
distribution leads to the detection of many small, artificial objects 
clustered around the true halo. Since these artefacts tend to be tran- 
sient, we perform an unbinding calculation using the SUBFIND al- 
gorithm (Springel et al. 2001; Dolag et al. 2008), and omit from our 
analyses any FoF haloes that do not have at least one self-bound 
substructure. This additional procedure successfully removes arti- 
ficial haloes from our sample. The version of SUBFIND used here 
is modified from the standard implementation, such that baryonic 
particles are also considered when identifying self-bound substruc- 
tures (for more details see Dolag et al. 2008). This provides an un- 
ambiguous definition of a galaxy within the simulations as the set 
of star particles bound to individual substructures. Thus, an indi- 
vidual halo may host more than one galaxy. 

When quoting halo masses, we refer to the total mass of the 
FoF halo in all components (i.e. gas, stars and cold dark matter). 
In some cases, however, it is appropriate to modify this definition 
slightly. When comparing results with simulations that follow only 
dark matter (e.g. Section 3), we consider only the cold dark mat- 
ter component of a FoF halo and boost its mass by a factor of 
f2m/(r2m — fib) to account for the baryonic component. Addi- 
tionally, when measuring baryon fractions it is common to specify 
fractions within the well-defined volume of a spherical overden- 
sity (SO) group (e.g. Lacey & Cole 1994); in this case we quote 
the mass within a sphere, centred on the local minimum of the 
gravitational potential, whose radius encloses a mean density (in 
all components) of 200pc, where rhoc is the critical density for a 
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flat Universe. We also quote the circular velocity, Vc, of haloes at 
this radius, r2oo. 



3 EVOLUTION OF THE HALO AND STELLAR MASS 
FUNCTIONS 

3.1 The halo mass function in different environments 

As pointed out by Frenk et al. (1988), and discussed more recently 
by Cole (1997) and Sheth & Tormen (2002), a fundamental differ- 
ence between regions of varying overdensity is the rate at which 
their dark matter halo populations evolve. This is reflected in the 
clustering bias of haloes as a function of mass (Mo & White 1996). 
As we shall see below, many environmental effects revealed by our 
simulations, for example, variations in the star formation rate den- 
sity, can be traced back to differences in the halo mass function in 
different environments. We therefore preface the discussion of star 
formation with a brief overview of the evolution of the halo mass 
function in each region, including resolution tests. 

In Fig. 2, we show the differential number density of haloes 
identified with the FoF algorithm, dn{M) /dM, multiplied by 
in order to reduce the dynamic range of the plot, and thus more 
clearly highlight the differences between regions. As noted in Sec- 
tion 2.3, the halo mass in this plot is taken to be Mpop = 
M°o|^f2m/(0,n — so as to enable the most direct compari- 
son possible with the results of dark matter only simulations. The 
regions simulated at intermediate resolution are shown as coloured 
curves, while the high-resolution realisation of the —2a region is 
shown as a grey curve. In the region of overlap, the agreement be- 
tween the high- and intermediate-resolution realisations is excel- 
lent, suggesting that the mass function in the latter has converged 
down to 1 X 10^" h^^ Mq (corresponding to ~ 200 particles). Ex- 
trapolating to the high-resolution simulation, we estimate that the 
halo mass function is reliable down to ~ 10® Mq. For subse- 
quent analyses of the baryonic properties of haloes, we therefore 
consider only FoF haloes with at least 200 particles. 

The GiMic mass functions are also compared with the global 
mass function fit of Reed et al. (2007), which extends to both 
higher and lower mass scales since it was derived from an ensem- 
ble of dark matter simulations covering a wide dynamic range in 
mass (including the Millennium Simulation; see also Jenkins et al. 
2001). As expected, the mass functions of the under- and overdense 
GiMiC regions bracket the global function, which is reasonably 
well matched by the mass functions of the and +lcr regions. We 
have checked that the halo mass functions of the regions al z = 1.5 
(when they were selected) are consistent with those of all possi- 
ble candidate spheres in the Millennium Simulation. As an exam- 
ple, the inset shows the multiplicity function of the intermediate- 
resolution Oct region, overplotted on the locus that encompasses 
the le**" and 84*'^ percentiles of the multiplicity functions of 100 
spheres, drawn from the Millennium Simulation, that quaUfied as 
candidates for the Ou region (for details see the Appendix). 

Although the different halo mass functions appear as scaled 
versions of one another at high redshift, a large difference devel- 
ops over time, as the amplitude and shape of the function evolves. 
By the present day, there are approximately 3 times more haloes of 
intermediate mass in the +2a region compared to the — 2cr; the dif- 
ference increases with mass, with large mass haloes found only in 
the +2a region. To determine the efficiency of halo formation per 
unit mass, we show in the right-hand panel of Fig. 2 the halo mass 
function now normaUsed by the total enclosed mass, rather than 



by the total enclosed volume. Normalised this way, the differences 
between the GiMic regions are smaller but there is still a large sys- 
tematic variation from region to region at all masses. Haloes form 
more efficiently in the high-density regions because they are the 
most dynamically advanced and the most massive haloes form only 
in the most overdense regions. These results agree with previous 
numerical and analytical studies of dark matter halo evolution (e.g. 
Frenk et al. 1988; Mo & White 1996; Sheth & Tormen 2002). 

3.2 The galaxy stellar mass function 

The galaxy stellar mass function as a function of redshift provides 
a useful check of the realism of our simulations. This function is 
shown in Fig. 3 and compared with observational data at « = 2 
and z = 0. The stellar mass functions in each of the five GiMiC 
regions at intermediate resolution are shown as coloured curves. 
From these, we obtain an estimate of the stellar mass function for 
the entire Millennium Simulation volume (black cun'e), by a suit- 
ably weighted average; details of this procedure are given in the 
Appendix. 

As in Fig. 2, we also plot results from the high-resolution re- 
alisation of the —2a region (grey curve) as a numerical test. This 
shows that at z = 2, the intermediate-resolution stellar mass func- 
tion has approximately converged for M, ^ 10^ Mq. Be- 
low this value, there are roughly twice as many galaxies in the 
intermediate-resolution simulation. This excess almost certainly re- 
flects a reduction in the efficiency of supernovae feedback in poorly 
resolved galaxies. At z = 0, the convergence properties are similar: 
below Mi, ~ l(f M©, the intermediate-resolution simulation 
overpredicts the number of galaxies by about a factor of 2, and 
above this mass, it underpredicts it by a similar factor. 

The weighted average stellar mass function closely tracks that 
of the Oa region over the mass range for which haloes are present in 
this region. At the high mass end, only the overdense regions con- 
tribute to the stellar mass function. This highlights the need to fol- 
low a wide range of large-scale environments when attempting to 
compare simulations with large galaxy redshift surveys. For com- 
parison, we overplot the stellar mass functions derived from the 
FORS and GOODS deep fields at 2 = 2 (Drory et al. 2005) and 
from SDSS data at 2 = (Li & White 2009). At 2 = 2, the data 
span a factor of 100 in stellar mass. Over this limited range, the 
model stellar mass function has a similar shape to the data but the 
amplitude of the weighted model function is approximately 0.5 dex 
higher, which is a slightly less than the difference between this 
model function and that of the —2a region. 

At 2 = 0, the data span nearly 4 orders of magnitude in stel- 
lar mass. The simulations produce too many small galaxies, with 

< Vf Mq, and too few in the mass range 10® - 10^° Mq, 
but inspection of the cumulative function shows that they produce 
about the right number density at M* > 10® M© and then at 
M* > 2.5 X 10^° Mq. At larger masses, the two overdense GiMiC 
regions produce a population of massive galaxies that is not seen in 
the data. The overproduction of faint galaxies and the "dip" at in- 
termediate masses result from the simple wind model that we have 
adopted. The excess at large masses in the overdense regions results 
from our neglect of heating processes that can quench cooling flows 
in clusters. In semi-analytic models of galaxy formation (Bower 
et al. 2006; Croton et al. 2006; Somerville et al. 2008) and cosmo- 
logical hydrodynamical simulations (Sijacki et al. 2008; Di Matteo 
et al. 2008; Booth & Schaye 2009) this is achieved by invoking 
feedback from AGN, which is not modelled in our simulations. 

The sensitivity of the stellar mass function to variations in the 
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Figure 2. Le/f/ Differential number density of haloes as a function of mass at z = 0, normalised by volume, and multiplied by M"^ in order to reduce the 
dynamic range of the plot. Coloured lines represent the intermediate-resolution GiMIC regions, and the grey line represents the —2a region at high resolution. 
The global halo mass function of Reed et al. (2007, black) is bracketed by the extreme GiMIC regions. Inset is the halo mass function at z = 1.5, plotted for 
the Oct region. The green shaded area delineates the 16^^ and 84"^ percentiles of the mass functions of 100 candidate Oct spheres drawn from the Millennium 
Simulation. Right: As left, but the differential number density is normalised by the total enclosed mass of each region, rather than by its volume. The difference 
between the regions is smaller when the regions are normalised by mass, but there is still a systematic variation from region to region at all masses. The most 
massive haloes form only in the most overdense, and hence most dynamically advanced, regions. 




7.0 8.0 9.0 10.0 11.0 12.0 7.0 8.0 9.0 10.0 11.0 12.0 

log,o M. [h-' Mg] log,„ M. [h"' M^] 

Figure 3. The stellar mass function of galaxies at z = 2 (left) and z = (right). Results are shown for all five intermediate-resolution simulations {coloured 
curves), and their weighted average {black curve), which represents an estimate of the function over the entire Millennium Simulation volume. The stellar 
mass function of the — 2ct region at high resolution is also shown {grey cun'e) to illustrate the degree of convergence. Observational data from the FORS and 
GOODS deep fields (Drory et al. 2005) are overplotted at z = 2, and from SDSS DR7 (Li & White 2009) at z = 0. The simulations are in rough agreement 
with the data at z = 2; at z = the overall number of galaxies more massive than ~ 10^ Mq is coiTectly reproduced by the simulations, but the number 
density of intermediate mass galaxies is too low, whilst that of the most massive galaxies is too high. 



physical assumptions and parameters in the simulations will be ex- 
plored in the OWLS project (Haas et al., in prep). For our pur- 
poses, it is sufficient to note that there is reasonably good agree- 
ment between the simulations and the data at high redshift while, 
at low redshift, the simulations produce roughly the right number 
of galaxies, albeit with an incorrect distribution of stellar masses. 



4 THE EVOLUTION OF COSMIC STAR FORMATION 

The cosmic star formation rate density (SFRD) is a key tracer of 
structure formation in the universe. After more than a decade of ob- 
servational studies, it appears to be fairly well constrained at least 
out to intermediate redshifts {z < 3), given assumptions about the 
initial mass function and reddening (Hopkins 2007). Taking into 
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Figure 4. Left: Star formation rate density, pi,(z), as a function of redshift in regions of different overdensity {colours), at intermediate {solid lines) and high 
resolution {dashed lines). We have not run the +2(t simulation at high resolution. The vertical dotted lines denote the epochs of reionisation of hydrogen 
(Hi) and helium (Hell) respectively. The variation in amongst the large-scale environments is pronounced; the most extreme regions differ by an order of 
magnitude at all epochs. Right: The evolution of the total SFR per unit mass, in regions of different large-scale overdensity. The persistent offsets demonstrate 
that the increase of pi, with large-scale overdensity does not arise purely because overdense regions enclose a greater mass. 



account the mass fraction recycled during stellar evolution, an in- 
tegral over p*(z) gives the present day cosmic stellar mass density 
(e.g. Cole et al. 2001; Rudnick et al. 2003; Eke et al. 2005; Li & 
White 2009). These are all quantities that are calculated in our sim- 
ulations. 

In this section we investigate the star formation history in our 
simulations and how it varies from one GiMIC region to another. 
We discuss the numerical convergence of our results and the na- 
ture of the dominant contributors to the star formation activity at 
high-redshift. The averaging technique described in the Appendix 
is applied to construct an estimate of the global SFRD that can 
be compared both to observations and to the results of the semi- 
analytic calculations applied to the Millennium Simulation. 

4.1 Large-scale environmental modulation 

For each region we consider all gas particles within the sampling 
sphere discussed in the Appendix, and use the volume of the sphere 
to compute the star formation rate per unit volume, [z), shown in 
Fig. 4. Comparison of the intermediate- and high-resolution simu- 
lations shows that our results have not converged at 2 > 6; the 
high-resolution realisations exhibit significantly higher overall star 
formation rates than their intermediate-resolution counteiparts. At 
later times, however, convergence of (z) is achieved to within 
~ 50 per cent, and inspection of the curves for the —2a region re- 
alisations demonstrates that this conclusion extends to the present 
epoch. 

The evolution of p* (z) may depend on resolution because the 
simulation fails to form small mass haloes, and/or because the star 
formation rate in the haloes that do form depends on resolution. 
The particle mass of the high-resolution realisations was chosen so 
that the Jeans mass in the photoionised IGM is marginally resolved. 
Hence these simulations should be able to follow the formation of 
all haloes in which the gas cools by atomic line cooling after reioni- 
sation. We show later that the star formation rate for resolved haloes 



is very similar at the two resolutions. We are therefore reasonably 
confident that our estimate of /6*(z) is close to numerically con- 
verged in our high-resolution simulations for z < «reion ~ 9, but 
clearly this is not the case in the intermediate-resolution realisa- 
tions for z > 6. The actual value of the star formation rate is, of 
course, sensitive to our subgrid modelling, in particular to the im- 
plementation of galactic winds. Note that the lack of high-redshift 
agreement in the value of between the two resolutions has only 
a small effect on the total amount of stars formed and on the dis- 
tribution of stellar ages at low-redshift, because the duration of the 
high redshift period before the values converge is small compared 
to the age of today's universe. 

The main features of the star formation rate density in our 
simulations are: 

• p* increases with decreasing redshift, peaks at 2: ~ 2 — 3, and 
then drops rapidly by a factor ^ 6 to 2; = 0; 

• the shape of p* (2) is similar for all regions; 

• the amplitude of Pt{z) varies by an order of magnitude be- 
tween the most extreme regions, at all epochs. 

The amplitude increases monotonically from the —2a to the +2a 
region at all epochs, but the difference between the +la and the 
+2(7 regions is far smaller than the difference between the —2a 
and the — Icr regions. The large variations in the amplitude of p* (2) 
between the different regions is not simply caused by the greater 
mass contained in the higher density volumes. As shown in the right 
hand panel of Fig. 4, higher a regions also have a higher specific 
star formation rate, Afstar/A/totai (where Aftotai is the total mass 
enclosed in each sphere). The SFR per unit mass still varies by 
~ 0.5 dex between the regions. We show later in Section 5.2 that 
this is due to the different halo mass functions in regions of different 
overdensity, and the variation in star formation efficiency with halo 
mass. 

A notable feature of the SFR in the high-resolution simula- 
tions (dashed lines) is the kink in p* at 2 ~ 9, the epoch when 
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Figure 5. Left: Star formation rate density, Pi,{z), in the L = 500 Mpc Millennium Simulation volume (black solid line), derived from a weighted 
average of the intermediate-resolution GiMIC simulations. Coloured lines show the instantaneous contribution to the total from dwarf galaxy haloes (blue), 
regular galaxy haloes (green) and group/cluster haloes (red). The evolution of predicted by the semi-analytic calculations of Bower et al. (2006, black 
dashed line) and De Lucia et al. (2006, black dot-dashed line) are shown for comparison. Symbols with error bars are observational estimates from the 
compilation of Hopkins (2007), converted to the Chabrier IMF. Right: The "archaeological" description of the star formation rate density. Haloes are assigned 
to bins according to their mass at 2: = and the coloured curves show the value of (z) due to their progenitors. 



the IGM is reionised and heated. The rapid increase in prior 
to this epoch is slowed by reionisation in all regions. The onset 
of reionisation has two main effects on the evolution of gas: (i) 
the photoheating boosts and maintains the gas temperature at a 
near-uniform level of ~ 10* K, and (ii) the photoionisation sup- 
presses cooling due to line emission from hydrogen and helium (Ef- 
stathiou 1992) as well as heavy elements (Wiersma et al. 2009a). 
Small haloes, whose virial temperature is below the temperature 
of the photoionised IGM, cannot hold on to their photoheated gas 
at 2 < Zroion, and star formation within them ceases (Pawlik & 
Schaye 2009). Both Hoeft et al. (2006) and Okamoto et al. (2008) 
used numerical simulations to construct a simple analytic model for 
how the characteristic mass, Mc(z), below which haloes lose most 
of their baryons, depends on the shape of the temperature-density 
relation for the IGM. This characteristic mass is not well resolved 
in the intermediate-resolution simulations (solid lines) for the ion- 
ising background adopted here (Haardt & Madau 2001), which ex- 
plains the absence of the kink in that simulation. In accord with 
the star formation history shown in the left-hand panel of Fig. 5, 
this also indicates that star formation at this epoch is dominated 
by low-mass haloes (Tvir < 10* K); otherwise reionisation would 
have had little effect on the global star formation rate density. 

The maximum in the star formation rate density at z ~ 2 — 3 
occurs significantly later than in the simulations of Springel & 
Hemquist (2003b), which have comparable numerical resolution 
to ours (and only slightly different cosmological parameters; they 
adopted = 0.3, VIa = 0.7, = 0.04, = 0.9, h = 0.7). 
We believe that this is only partly due to our inclusion of metal 
line cooling that affects the SFR particularly at lower redshift, as 
discussed and illustrated in the analytical model of Hernquist & 
Springel (2003). The largest simulations of Oppenheimer & Dave 
(2006) cover volumes of (32 Mpc)'^ and have 256'^ particles, 
which results in a factor of 2 poorer mass resolution than our 
intermediate-resolution simulations, and a factor ~ 15 poorer res- 



olution than the high-resolution GiMIC runs. At this level of reso- 
lution, our simulations would certainly be far from converged. The 
SFR in most of their models also peaks at a higher redshift than 
ours. The simulations of both Springel & Hemquist (2003b) and 
Oppenheimer & Dave (2006) use the multiphase ISM implementa- 
tion of Springel & Hernquist (2003a) for gas in galaxies, as well as 
a different prescription for the generation of galactic winds. The de- 
tails of the ISM and wind implementations are probably the main 
cause of the formation history differences between these simula- 
tions and those presented here. 



4.2 The Millennium Simulation star formation rate density: 
comparison with observations and semi-analytic models 

The star formation rates in the individual GiMIC regions can be 
combined in the manner discussed in the Appendix in order to es- 
timate the net star formation rate of the whole Millennium Simu- 
lation volume. In order to obtain results down to z = 0, we use 
the intermediate-resolution simulations; we consider only epochs 
z < 6 where our results are approximately converged. 

Our estimate of p* (z) for the Millennium Simulation volume 
is shown in Fig. 5. It increases with time to reach a broad plateau 
atRi O.25M0yr-iMpc-^ over the interval 3 > z > 1 with a 
maximum around z ~ 2, followed by a rapid decline by a factor of 
~ 6 to z = 0. The behaviour of the global p* (z) closely follows 
that of the mean density (Oa) region. 

The coloured curves in the left-hand panel decompose p* (z) 
into contributions from haloes binned according to their mass 
at that epoch; bins are chosen to correspond approximately to 
the masses of haloes that host dwarf galaxies ("dwarf haloes"; 
M < 5 X lO^^/i^^Mo), typical galaxies ("galaxy haloes"; 
5 X 10" ft-^ Mq < M < 5 X 10^^ Mo) and groups and 
clusters of galaxies ("group haloes"; A/ > 5 x lO'^^ Mq). As 
might be expected within the CDM hierarchical assembly model. 
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the peak of star formation within low mass haloes precedes that of 
more massive haloes. For our particular classification, these peaks 
are broad, occurring at z ~ 4 — 5 for dwarf haloes, z ~ 2 — 3 for 
galaxy haloes, and 2 ~ 1 — 2 for group haloes. 

The simulation outputs record the. formation time of each star 
particle and this allows us to relate present day systems to the prop- 
erties of the progenitors in which their stars formed. This "arche- 
ological" description of the star formation history is computed for 
haloes of a given mass at 2 = by summing the initial mass of their 
stellar particles (i.e. the mass at the time of formation, prior to mass 
loss due to stellar evolution), binned by their formation redshift, 
and is shown in the right-hand panel of Fig. 5. Note that the solid 
black curve, integrating over haloes of all masses, differs slightly 
from that of the left-hand plot at epochs other than z = Q, because 
the population of baryonic particles within the sphere used to de- 
fine the volume limited "sample" for these plots (see the Appendix 
for details of how this is selected) evolves slightly over time. This 
plot demonstrates that star formation is dominated at all redshifts 
by the progenitors of galaxies that today reside in the most massive 
haloes (i.e. groups and clusters). This is perfectly compatible with a 
hierarchical build-up, because the star formation in massive haloes 
became dominant only recently, as shown in the left-hand panel of 
the figure. 

To compare the GiMIC simulations to observations, we use 
the compilation of star formation rates given by Hopkins (2007). 
For consistency, we have adjusted the published values to account 
for the cosmological parameters and IMF adopted in the simula- 
tions; to convert from the IMF of Salpeter (1955) to that of Chabrier 
(2003), pi, is scaled down by a factor of 1.65. As noted in Section 
2.1, the mass loading in the galactic winds was chosen by requiring 
small test simulations to produce an approximate match to the am- 
plitude of the observed p-k{z). As a consequence, the broad agree- 
ment with the data is largely a check that the weighting scheme 
works well. 

The evolution of the star formation rate density shows a broad 
plateau of nearly constant star formation rate during 4 > 2 > 1, 
followed by a rapid drop to 2 = 0. The reduction in towards 
the present day is smaller in the simulation than in the data. This 
is mostly due to the high star formation rates in massive haloes in 
the simulations (red curve in the left-hand panel of Fig. 5), which 
is also the cause for the overproduction of bright galaxies seen in 
Fig. 3. Recent semi-analytic models (e.g. Bower et al. 2006; Croton 
et al. 2006; De Lucia et al. 2006; Font et al. 2008) and cosmological 
hydrodynamical simulations (Sijacki et al. 2008; Di Matteo et al. 
2008; Booth & Schaye 2009) include AGN radio-mode feedback 
to suppress the cooling of gas in such haloes and quench their star 
formation. As a result, p* in these models falls faster from 2 = 1 
to 2 = 0, in better agreement with the data (dashed and dot-dashed 
lines in Fig. 5). Kobayashi et al. (2007) obtained a qualitatively 
similar star formation history to that of GiMIC in hydrodynamic 
simulations with very strong feedback from hypemovae (HNe), but 
no AGN feedback. Like us, they found star formation at late times 
to be dominated by massive blue galaxies. 

The discrepancy in p* occurs at late times when the star for- 
mation rate in all objects is relatively low, and therefore has little 
effect on the global stellar mass density, SI* (in units of the critical 
density). The simulation produces a current stellar mass density 
(obtained by integrating p*(2:) and subtracting the mass recycled 
by stellar evolution) of f2* ~ 3.2 x 10"'^ at 2 = 0. Observational 
estimates of this quantity give fi* — (1.37 ± 0.4) x 10^** from the 
2dFGRS (Eke et al. 2005) and = (1.44 ± 0.4) x 10"^ from the 
SDSS (Li & White 2009), assuming the Kennicutt and Chabrier 
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Figure 6. The ratio of past average star formation rate to the current star for- 
mation rate at the redshift of observation, for three mass bins. The binning, 
matched to Fig. 7 of Bower et al. (2006), is done according to the stellar 
mass at the redshift of observation. The dotted line at p* /tn = P* sep- 
arates the regimes where i) on-going stai' formation dominates the cun'ent 
stellar mass (below), and ii) on-going star formation contributes a negligible 
fraction to the total stellar mass (above). 

IMFs respectively. These estimates have a systematic uncertainty 
of approximately a factor of 2 arising from the choice of IMF. 

4.3 Cosmic downsizing 

Observations indicate that star formation in the most massive galax- 
ies observed today essentially concluded at 2: ~ 1 (e.g. Drory et al. 
2003, 2005; Pozzetti et al. 2003; Fontana et al. 2004; Kodama et al. 
2004), and that their star formation rates were higher in the past 
than they are at present (e.g. Cowie et al. 1996; Juneau et al. 2005). 
These findings have been cited as a challenge to the ACDM cos- 
mogony, since a naive view of galaxy formation within a hierarchi- 
cal assembly framework would have star formation shifting from 
smaller to larger objects as the hierarchy builds up. This appar- 
ent conflict has spawned terms such as "downsizing" and "anti- 
hierarchical galaxy formation". Tentative evidence that these obser- 
vations are compatible with ACDM was first presented by Pearce 
et al. (2001), using hydrodynamic simulations that featured star for- 
mation but no feedback mechanisms. Bower et al. (2006) and De 
Lucia et al. (2006) convincingly demonstrated that the observations 
are entirely compatible with their ACDM semi-analytic models. 
We now explore whether our hydrodynamic simulations also ex- 
hibit downsizing. 

One measure of downsizing, first highlighted by Bower et al. 
(2006, their Fig. 7), is the ratio of the past average SFR of galaxies 
of a given mass, Mi,[z) /t-a^z), where in is the Hubble time, to 
their SFR at the epoch of observation M«(M,t, z). This is shown 
in a volume-averaged sense in Fig. 6, for which galaxies have been 
drawn from all five intermediate-resolution regions. Galaxies enter 
the regime above the horizontal dotted line that marks p*/tH = 
pi, when they have assembled most of their stellar mass and their 
present star formation is adding only a minor contribution. 

In common with results from the semi-analytic model of 
Bower et al. (2006), galaxies in GiMIC exhibit a segregation by stel- 
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Figure 7. Cumulative number density of galaxies as a function of stellar 
mass at 2 = 6 for the five GiMIC regions (solid coloured lines) at interme- 
diate resolution. The number density for the —2a region is also shown at 
2 = 3 (dark blue dashed line). The data plotted clearly illustrate the strong 
bias of massive stellar systems towards overdense regions. The most mas- 
sive galaxies, drawn from the +2cr region, have solar abundances (black 
points, right-hand y-axis) and are highly overabundant in a-process ele- 
ments (red points, right-handy-axis). 



lar mass, such that the least massive galaxies continue to grow until 
the present epoch, whilst more massive galaxies have essentially 
concluded their star formation at z > 1 — 2. Although this effect 
is weaker in our simulations than in the semi-analytic model and in 
the observational data, this is a significant result because our model 
does not include AGN radio-mode feedback. Rather than causing 
the downsizing effect directly, AGN feedback appears merely to ex- 
acerbate it, pushing massive galaxies further into the regime where 
the past average star formation rate dominates the current star for- 
mation rate. AGN feedback does seem to be required, however, to 
prevent the formation of the overly bright and blue galaxies that 
form in large haloes in our simulations but which are not observed 
in nature. The sort of downsizing discussed in this section results 
not from AGN feedback but rather from the complex interplay of 
gas cooling, star formation and feedback that develops as structure 
assembles hierarchically. 



Since the mean thickness of a sphere is 4r / 3, the comoving "depth" 
of each region is ~ 24 Mpc. The expansion rate at z = 6 is 
Ho = 930 /i km s"^ Mpc"\ yielding a mean velocity range of 
Aw ~ 3190 km s~^. This corresponds to Az ~ 0.07 and there- 
fore to a redshift "window" of z = 6 ± 0.035. Although this is 
much narrower than the Az = 1.5 estimated by Bouwens et al. 
(2008) for their survey at z = 7, the actual volume probed with 
the pencil-beam geometry of the HUDF at these epochs is similar 
to that of each GiMIC region. Our simulations indicate that p* (z) 
varies by up to an order of magnitude between regions of this vol- 
ume and so sample variance is a significant source of systematic 
uncertainty in measurements of at high redshift. This is further 
compounded by the dominant contribution to p* at these epochs of 
low mass, inefficiently star-forming galaxies that are below current 
detection limits. 

The high-resolution simulations produce a slowly varying star 
formation rate density from z = 6 to z = 9. If the escape fraction 
of ionising radiation from the small haloes that dominate p* at these 
epochs is high, say ~ (25 — 80) per cent as suggested by Wise & 
Cen (2009), then these galaxies will be the dominant contributors 
to the UV-background at high redshift, and thus may have been the 
main sources of the radiation that led to reionisation (Srbinovsky 
& Wyithe 2008). The faintness of each individual source would 
reconcile the apparent dearth of detected sources of UV-photons 
with the inferred ionisation state of the intergalactic medium at z ^ 
6 (Bolton & Haehnelt 2007). 

Although small galaxies dominate the star formation rate den- 
sity, the wide range of overdensities represented in the GiMIC re- 
gions results in massive stellar systems (i\f* ~ 10^^ M©) 
forming in the simulations as early as 2: = 6, as shown in Fig. 7. 
The dark matter haloes that host these massive galaxies are so 
strongly biased towards overdense regions that galaxies of com- 
parable mass do not form in the void region (—2a) until 2: ~ 3 
(dark blue dashed line). These massive galaxies are embedded in 
large, nearly spherical coronas ofhot(r~ 10^ K) gas of comovme 
radius ~ 0.3 h^^ Mpc, but the galaxies themselves are extremely 
compact (comoving radii of ~ 3 kpc). The stars in these galax- 
ies have near-solar metallicity and are highly overabundant in the 
a elements produced in type II SNe (such as oxygen), relative to 
those produced by type la SNe (such as iron). They are relatively 
old (ages ~ 0.2 Gyr) at this redshift. Such galaxies are reminiscent 
of the z > 5 massive and evolved galaxy candidates found in the 
GOODS fields by Wiklind et al. (2008). 



4.4 High-redshift star formation 

The observational data for 2: > 6 are uncertain and differ by up 
to a factor of ten between different studies; compare for example 
the results of Bunker et al. (2006) and Bouwens et al. (2008). Our 
high-resolution simulations favour a modest drop of ~ 50 per cent 
in p* from its peak at 2 ~ 2 to z = 6, with most stars forming in 
faint galaxies, which have star formation rates of < 0.1 Mq yr^^. 
Such star formation rates are typically an order of magnitude below 
the sensitivity limits of most current surveys. 

The sky coverage of each r = 18 Mpc GiMIC region is 
320 arcmin"^ at 2 = 6. Whilst this is comparable to the sky cover- 
age of the GOODS fields (Dickinson et al. 2003), it is much greater 
than the 5-20 arcmin^ of the Hubble Ultra Deep Field (HUDF, 
Beckwith et al. 2006) in various passbands, which is an important 
source for the detection and characterisation of galaxies at 2 > 6. 



5 THE STAR FORMATION PROPERTIES OF HALOES 

We now turn to an exploration of the physical basis for the environ- 
mental variation in the cosmic star formation rate density described 
in Section 4. We show that large-scale environment has no impact 
upon the evolution of the star formation rate in individual haloes. 
Instead, the strongly differing amplitudes of p* in the GiMIC re- 
gions can be traced back to their different mass functions. 

We begin by considering the star formation rate in haloes of 
different circular velocity and how this depends on the "subgrid" 
physics assumed in the simulations. We compare the star forma- 
tion properties of the simulated haloes with analytic models for the 
regulation of star formation by gas cooling, and show that winds 
effectively regulate star formation up to a threshold value of the 
halo circular velocity. In common with conclusions drawn in Sec- 
tion 4, from comparison with a semi-analytic model that includes 
AGN feedback, we infer that efficient feedback in massive, hydro- 
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static haloes is necessary to regulate the growth of the most massive 
galaxies and to reconcile the cosmic star formation rate density at 
low redshift in the simulations with observations. 



5.1 The specific star formation rate in haloes 

To characterise the star formation efficiency of haloes, we define 
the halo specific star formation rate, sSFR=A/*/A/2()0, where M200 
is the total mass (baryons plus dark matter) contained in a sphere 
of radius r ^oo- For clarity, we shall retain the "halo" prefix to avoid 
confusion with the more common definition of specific star forma- 
tion rate that applies strictly to galaxies, Mi,/Mi,. 

In common with previous sections, we start by discussing the 
degree of numerical convergence in the halo sSFR achieved by our 
simulations. The mean halo sSFR as a function of circular veloc- 
ity, Vc = {GM2oo/r2ooY^^, and M200 for haloes identified in the 
five intermediate-resolution regions is shown at various redshifts in 
Fig. 8. The small box and whisker elements denote one and two 
standard errors on the mean. Since the relation for each region is 
virtually identical at each redshift, convergence can be assessed by 
comparison with the relation from the high-resolution realisation 
of the —2a region for all epochs (dark blue solid and grey dotted 
curves). 

For Vc > 250 km s~^, the halo sSFR in the —2a re- 
gion is very similar for the two resolutions; for 100 km < 
"Uc ^; 250 km s~^, the high-resolution simulation is only slightly 
higher than the intermediate-resolution simulation. However, be- 
low 100 km s~^, the intermediate-resolution simulation clearly 
overpredicts the halo sSFR (except at the highest redshift) and be- 
comes approximately independent of Vc. Extrapolating from this 
behaviour to haloes resolved with the same number of particles in 
the high-resolution simulations, we conclude that these should give 
reliable results for Vc 50 km s^^. 

The halo sSFR-iic relations for all GiMIC regions are virtually 
identical at every redshift. This similarity may be unexpected, given 
the large variation in the between the regions; we defer to Sec- 
tion 5.2 a detailed explanation of how these results are reconciled, 
in order to focus initially on the origin of the halo sSFR-Vc rela- 
tion. The relation varies strongly with redshift, but drops rapidly 
below a characteristic circular velocity that is essentially indepen- 
dent of redshift. This drop occurs at 250 kin s^^, in the regime 
where the high- and intermediate-resolution simulations agree well, 
and is caused by the action of galactic winds. This scale is signif- 
icantly lower than the speed with which galactic winds are actu- 
ally launched {v^ = 600 km s^^, indicated by an arrow) because 
the effective velocity of the winds is suppressed by drag forces im- 
parted by the dense gas of the ISM. Dalla Vecchia & Schaye (2008) 
demonstrated that the winds become ineffective if the kick velocity 
falls below a critical value that increases with the ISM pressure and 
hence with halo circular velocity (or mass). Previous studies such 
as Springel & Hemquist (2003b) did not find this because the wind 
particles were temporarily decoupled from the hydrodynamics and 
hence could escape freely without experiencing any drag from the 
ISM. The feature becomes more pronounced at later times when 
the halo sSFR drops faster with time for lower masses. 

Haloes of virial mass M200 ~ lO^/i"^M0 and circu- 
lar velocity ~ 50 km s"'^ (that are well resolved in the —2a 
high-resolution simulation) have low star formation rates, M* < 
0.05 Mq yr^^, at redshifts 2; ~ 6. Such objects are ubiquitous at 
that time and they dominate the total star formation rate. The halo 
sSFR drops rapidly at lower masses as haloes lose their baryons due 
to a combination of stellar feedback and heating of their gas by the 



UV-background. The GADGET2 simulations of Springel & Hem- 
quist (2003b) yield quantitatively similar results, as do the AMR 
simulations of idealised dwarfs described by Wise & Cen (2009). 
In the latter, blow-out of baryons prevents continued star formation 
at masses below 10* Mq . In our simulations, the halo sSFR in these 
small haloes drops rapidly by more than two orders of magnitude 
towards z = 0. 

Above a circular velocity of ~ 250 km s~^, where haloes are 
less susceptible to galactic winds in our model, the halo sSFR is 

^ 10^^ Gyr^^, almost independently of mass above 2; ~ 3. The 
level of this plateau drops by a factor ~ 5 to z = 1. The simulations 
of Springel & Hemquist (2003b) yield a sSFR for a halo of mass 
10^^ h'^ Mq of 10"^-^ Gyr"^ at 2 = 3, a factor of 3 below our 
value. They also show the halo sSFR at z = 1.5, and the drop from 
a = 3 is a factor of 5, similar to the drop in our simulations over the 
interval z = 1 — 3. Whereas the halo sSFR is nearly constant with 
mass for massive objects in our simulations, it increases rapidly 
with halo mass in those of Springel & Hernquist (2003b). 

The halo sSFR plummets by a factor ~ 50 for a 10^^ Ir'^ Mq 
halo between z = 6 and z = 0, and by an even greater factor 
at lower masses. The cause of this strong redshift evolution can be 
understood from simple physical arguments. White & Frenk (1991) 
present an analytic model in which the halo sSFR of a massive ob- 
ject of given circular velocity drops with redshift as star formation 
becomes limited by the gas cooling time rather than by its free-fall 
time. Defining the cooling rate, A(m), via 

pu = -AriH, (3) 

where nn = Xp/mn is the hydrogen number density for gas at 
density p, of which a mass fraction X is hydrogen, and u is the 
thermal energy per unit mass, this cooling time is then 

u u u 



Tc(p) 



(4) 



M ^n^/p Ap 

The cooling function depends on temperature, T. The virial tern 
perature of a halo is related to its circular velocity by: 

T,i^= ^2^3.6 X 10^ K 



(0.59) ( 



100 kms 



3 kB 

where is the mean molecular weight of the gas. 

At high redshift, star formation is limited by the gas con- 
sumption timescale since Tc becomes very small. So long as in- 
verse Compton cooling is unimportant, which is true at low red- 
shifts, the cooUng time (and hence the halo specific star formation 
rate) of haloes with virial temperature ~ 10^ K (corresponding 
to Vc = 530 km s^^), scales oc A^^p^^ which is approximately 



p-' = (l + . 



if A does not evolve strongly. 



The evolution of the normalisation of the halo sSFR-Vc rela- 
tion according to this model was explored analytically by Hemquist 
& Springel (2003). Starting from pseudo-isothermal initial gas pro- 
files, their model yields specific star formation rates for haloes of 
Tvir — 10^ K that closely trace those in the simulations of Springel 
& Hemquist (2003b) at intermediate redshifts (7 > z > 2). This 
model is compared to the high-resolution GiMiC simulations in 
Fig. 9 where we have computed the cooling rate taking into ac- 
count the effect of an evolving ionising background as given by 
Haardt & Madau (2001) and described by Wiersma et al. (2009a). 
Note that this plot does not represent the evolution of the halo sSFR 
of a given object, but rather that of objects of a given circular ve- 
locity at different redshifts. 

The effect of metallicity on the sSFR is pronounced in the 
model, but interestingly we find that there is, in fact, very little 
metallicity evolution for the hot gas in the simulations, Z ~ 0. 1 Z© 
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Figure 8. The "halo specific star formation rate", M4./M2oo> as a function of circular velocity (v^, bottom axis) or halo mass (A/2oo> top axis) for various 
redshifts. Different colours are for each intermediate-resolution GiMIC region as indicated. Dotted grey lines are the results of the high-resolution —2a 
simulations. Box and whisker elements indicate one and two standard errors on the mean for the Otr region. The halo sSFR— iic relations of the five GiMIC 
regions are essentially indistinguishable and show a break at Vc ~ 250 km s~ ^ , independently of redshift. The redshift evolution of the halo sSFR either side 
of this break reflects the differing mechanisms acting to regulate star formation at different halo mass scales. The arrow at 600 km s^^ indicates the launch 
speed of the winds. 



from 2 = 9 to the present. For this value of the metallicity, the 
model works very well for massive haloes, indicating that the cool- 
ing time controls star formation in these objects in the simulations, 
as expected from the White & Frenk (1991) calculation. (Note that 
the inclusion of an efficient feedback mechanism in massive haloes, 
such as energy injected by AGN, could alter the sSFR). The sim- 
ulation results and the Hernquist & Springel (2003) model diverge 
at high redshift because these authors chose not to include the gas 
consumption timescale in this calculation; however, as White & 
Frenk (1991) show, taking account of this constraint would reduce 
the star formation rate and improve the match at high redshift. 

Efficient feedback depresses the halo sSFR in low-mass ob- 
jects far below the rate expected from the cooling time argument 
alone. The accelerated drop-off (relative to the trend for massive 
haloes) at z < 3 is due to feedback from galactic winds which de- 
pletes the baryon fraction in small haloes and increases the cooling 
time of the remaining gas, quenching star formation. Winds there- 



fore introduce a second characteristic scale in galaxy formation in 
our model, that at which star formation is quenched due to the ejec- 
tion ofbaryons. For our chosen wind speed of = 600 km s~^, 
this scale is around Vc ~ 250 km s^^ and depends on the details of 
the wind implementation. This value is close to that advocated by 
White & Frenk (1991) as the scale separating massive haloes with 
long cooling times from smaller haloes within which gas cools ef- 
ficiently. 



5.1.1 Comparison with semi-analytic models 

Although we defer a detailed comparison of the GiMIC and semi- 
analytic treatments of galaxy formation in the Millennium Simula- 
tion to a later paper, it is instructive to carry out a limited compar- 
ison here. We focus on the halo sSFR, M*/M2oo, and compare, in 
Fig. 10, the results obtained by Bower et al. (2006) and De Lucia 
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Figure 9. The mean halo specific star formation rate, Af*/M2oo, as a func- 
tion of redshift for haloes of circular velocity t^c < 300 lim {solid 
red line) and = 50 km (solid blue line), drawn from the high- 
resolution — 2(7 GiMIC simulation. Results are compared with the ana- 
lytic model of Hemquist & Springel (2003, labelled HS03) for the massive 
haloes, assuming three different hot gas metallicities (primordial, dotted 
red; 0.1 Zq, dashed red; 1.0 Zq dot-dashed red). The metallicity of the 
hot gas in the simulations is nearly constant at ^ 0.1 Zq, and the model 
predicts the result for the massive halo in the simulation very well. At lower 
redshifts the drop in star formation follows that expected from the decrease 
in density. The analytic model is a poor description for low mass haloes 
(blue), since it neglects the action of winds; the predicted halo sSFR is or- 
ders of magnitude higher than the simulated values, and so is not shown on 
the plot. 



et al. (2006) with those from GiMIC. Note that the GiMIC sim- 
ulation extends to lower halo masses because of its higher mass 
resolution. 

For halo masses M200 > 10^^ h^^ Mq, the two semi-analytic 
models give very similar results, particularly for 2; < 3. At smaller 
masses, however, there are substantial differences between the 
semi-analytic models and also between these and the simulations. 
These differences are largely due to different treatments of feed- 
back. The de Lucia et al. model has substantially higher sSFR than 
the Bower et al. model reflecting the weaker feedback assumed 
in the former; the feedback adopted in GiMIC is intermediate be- 
tween the two semi-analytic prescriptions. For halo masses around 
3 X 10^^ Mq, the GiMIC results agree well with the semi- 
analytic models. However, at higher masses GiMIC produces much 
higher star formation rates. The reason for this discrepancy is sim- 
ply that the semi-analytic models include feedback from the AGN 
radio-mode. As explained in their respective papers (see also Ben- 
son et al. 2003; Croton et al. 2006), this "radio mode feedback" is 
required in order to prevent the formation of overly massive galax- 
ies in large haloes. 

Finally, dotted and dashed lines in the bottom-right panel of 
Fig. 10 compare the Bower et al. (2006) halo sSFRs in the +2(t 
and —2a GiMIC regions. As in the GiMIC simulations (see Fig. 8), 
the semi-analytic treatment finds little, if any, dependence of the 
halo sSFR on the large-scale environment. 
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Figure 10. The halo specific star formation rate. A/* /AI20Q, as a function of 
halo mass, M200, (as in Fig. 8) for all GiMIC regions combined, compared 
to the semi-analytic models of Bower et al. (2006, red lines) and De Lucia 
et al. (2006, blue lines). At z = 0, results of the semi-analytic model are 
shown, limited to the — 2cr (red dotted curve.s) and +2cr (dashed curve.s) 
regions, rather than the entire Millennium Simulation volume; no difference 
is apparent, as in the hydrodynamic simulations. 

5.2 The multiplicity function of star formation and its 
environmental dependence 

We now return to an explanation of the physical reason for the 
strong dependence of the cosmic star formation rate density, p*, 
on large-scale environment presented in Section 4.1. As shown in 
Fig. 4, p* varies by approximately an order of magnitude amongst 
the GiMIC regions, which may seem surprising given the similarity 
of the sSFR-«c relation in each region. 

The halo star formation rate density may be written in terms 
of the specific star formation rate per halo as: 



M21 



dN{M2oo, 



din M200 
g{M2oo, z)dlnM2oo 



M200 



dlnM2. 



(6) 



where dN{M2oo ,z)/dln M200 is the mass function of dark matter 
haloes. The function 



g{M2oo, z) = M200 



dN{M2oo,z) 



(7) 



din M200 M200 
is the multiplicity function of star formation, which describes the 
star formation rate density in haloes of a given mass. As we showed 
in Section 3, the halo mass function is a strong function of environ- 
ment, with the largest haloes forming only in the densest regions. 
However, as we showed in Section 5.1 the star formation properties 
of a halo at a given epoch, including its specific star formation rate, 
i\f* /A/200, depend on its mass, not on its environment (Fig. 8), i.e. 
at any given epoch Af* /Af200 is essentially independent of large- 
scale environment. The dependence of on large-scale environ- 
ment therefore reflects the dependence of the halo mass function 
on large-scale environment. 
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Figure 11. The multiplicity function, g(M2oo, z), of star formation (shaded) as a function of region (rows) andredshift (columns). This function is the product 
of the halo sSFR and the halo multiplicity function, both of which are functions of redshift and halo mass. For reference, the halo multiplicity function of each 
region is reproduced from Fig. 2 {black histogram, right-hand axis), as is the multiplicity function derived from the universal mass function fit of Reed et al. 
(2007, dotted curve, right-hand axis). 



The multiplicity function of star formation, g{M2oo, z), is 
shown in Fig. 1 1 for each GiMIC region at three different redshifts; 
the star formation rate density of each region is proportional to the 
shaded area. As highlighted by Springel & Hernquist (2003b), the 
history of p* {z) is shaped by the build-up of matter in haloes, as de- 



scribed by dN/d In Af2oo, and the evolution of the specific star for- 
mation rate, M*/M2oo. On the one hand, advancing structure for- 
mation shifts an ever increasing fraction of bound mass into more 
massive haloes, thereby pushing more and more gas above the ef- 
ficiency thresholds imposed by the photoionising background and 
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galactic winds. Conversely, the amplitude of Nh/M'zoo is greatest 
at early times and drops steadily over time. The combination of 
these effects yields a broad plateau in at intermediate redshifts 
(6 ^ z ^ 2), from which it drops off towards lower and higher 
redshifts. 

Massive haloes are underrepresented in low density or "void" 
regions and consequently star formation within voids is dominated 
by low-mass haloes at high redshift (i.e. z = 6), with the multi- 
plicity function of star formation peaking at ~ l(f Mq. Only 
overdense regions have sufficient numbers of massive haloes at 
these early times for massive galaxies to contribute significantly 
to pi,; in these regions a second peak develops at ~ 10^^ M©. 
By 2 = 1, massive haloes begin to be found even in voids and 
so the first peak in g moves to haloes of mass ~ 10^^ h^^ Mq in 
all regions with a second one developing around ~ lO'^^ Mq 
in the high density regions. At the present time, star formation is 
dominated by galactic haloes (M200 ~ 10^^ Mq) in the void 
region, by groups (A/2o() ^ 10^^~^* Mq) in the intermedi- 
ate density regions, and by clusters (M200 10^ ' h^^ Mq) in the 
highest density region. It is the absence of massive, efficient star- 
forming haloes in voids that so severely inhibits their overall star 
formation rate density. This explains the much greater difference 
in between the underdense —2a and — la regions, compared to 
that between the overdense +la and +2a ones. 

Our simulations overestimate the specific star formation rates 
in massive galaxies and haloes; hence the contribution of clusters is 
likely to be exaggerated. Even so, the multiplicity function plot il- 
lustrates the dynamic range problem faced by simulations of galaxy 
formation: is dominated by 10^ Mq haloes at z = 6 but (in 
our simulations at least) by 10^® h^^ Mq haloes at z = 0. As an 
example, the Millenniimi Simulation does not resolve haloes below 
2 X 10^° h-^ Mq and hence misses the objects that (in our simu- 
lations) are the dominant contributors to star formation at z = 6. 
On the other hand, periodic simulations of, for example, side length 
L = 32 Mpc enclose a total mass of ~ 10^® Mq, and so 
cannot contain even a single cluster of the type that dominates p^ 
in the +2a region. However, the fact that M*/M2oo turns out to 
depend mostly on halo mass and not on overdensity makes it possi- 
ble, in principle, to obtain accurate values for pi, and its dependence 
on environment. This may be achieved by combining accurate esti- 
mates of M*/M2oo as a function of halo mass with the dependence 
of the halo mass function on environment. 

5.3 Halo baryon fractions 

Stars comprise a relatively small fraction of the total baryonic con- 
tent of the Universe. In this subsection, we investigate how the to- 
tal baryonic content of haloes varies with halo circular velocity and 
with time. This allows us not only to generate an inventory of the 
baryons in our simulations but also to understand further the pro- 
cesses responsible for the star formation properties discussed in the 
previous subsection. As was the case for the sSFR, the baryon frac- 
tions of haloes are essentially independent of environment and so, 
when necessary, we combine haloes from all 5 GiMIC regions to 
improve the statistics. As before, we begin this discussion with a 
convergence analysis. 

The baryon fraction in tmits of the mean, /b = 
(Mb /M2oo)/(f^b /f^m), is plotted as a function of circular velocity 
(and halo mass, top axis) at several redshifts in Fig. 12. To investi- 
gate convergence, we use the intermediate- and high-resolution ver- 
sions of the —2(7 region. The left-hand plot shows that approximate 
convergence is achieved at a circular velocity of Vc ~ 100 km s~^. 



the same value at which, as shown in the preceding subsection, ap- 
proximate convergence in sSFR is also achieved. As before, ex- 
trapolating from the behaviour seen in Fig. 12, we assimie that 
baryon fractions in the high-resolution simulation are reliable down 
to Vc ^ 50 km s^^ (that corresponds to the resolution limit scal- 
ing with the number of particles in the intermediate- and high- 
resolution simulations.) The rapid upturn in the baryon fraction at 
Vc < 30 km is thus Ukely to be an artefact of the limited res- 
olution. 

The right-hand plot of Fig. 12 shows how the baryon content 
of a halo is made up from the following components: 

• stars 

• gas with nn > 0.1 cm~^ (the value above which we impose 
an equation of state to mimic the multiphase ISM - see §2.2.) 

• cold gas 

• warm gas 

• hot gas. 

The cold and warm gas components are separated at T = 
10^ K, whilst the warm and hot gas components are separated at 
T = 10^'^ K. These limits are chosen so as roughly to distinguish 
the warm-hot medium from hotter X-ray emitting gas. The frac- 
tions of the different components are plotted cumulatively in the 
order listed above, so that the curve showing the hot component is 
equivalent to the total baryon fraction. For Vc < 160 km s~^, we 
consider only the (well resolved) haloes from the high-resolution 
—2a simulation; for larger values of Vc, we draw haloes from all 
five intermediate-resolution runs to obtain a large sample of well 
resolved haloes. 

Massive haloes with circular velocity comparable to the 
launch speed of winds or greater hold on to their baryons and, as 
shown in Fig. 12, have /b ~ 0.9, as expected in the non-radiative 
regime (e.g. Eke et al. 1998; Grain et al. 2007). The baryon frac- 
tion drops rapidly for haloes with Vc < 200 km s~^. These haloes 
have a virial temperature that coincides with the peak in the coohng 
function. Efficient gas cooling leads to high rates of star formation, 
but the associated galactic winds are powerful enough to drive the 
gas out of the potential wells. As a result these haloes have, in fact, 
lower specific star formation rates (see Fig. 8) and baryon fractions 
than more massive haloes. 

Whereas the total baryon fraction in haloes with Vc > 
300 km is nearly independent of redshift, the fraction in the 
star-forming "equation of state" gas phase drops dramatically from 
35 per cent at « = 6 to a negligible fraction at 2 = 0, when ~ 75 
per cent is in the form of hot gas. The overall baryon fraction of less 
massive haloes, Vc ~ 100 km s^^, drops significantly with red- 
shift. Their stellar fractions increase with time and become roughly 
constant for 2 < 3, at around 5 per cent of the cosmic mean. Be- 
cause of their short gas cooling time, most of the gas in these low 
mass haloes is initially in the cold phase. As star formation pro- 
ceeds, winds typically eject over half of the baryons that are par- 
tially converted into warm and hot gas. Other mechanisms such as 
ram-pressure stripping of gas and tidal stripping of stars can con- 
tribute to the removal of baryons from these haloes. 



6 SUMMARY AND CONCLUSIONS 

We have presented the first results of a programme of hydrody- 
namical simulations of the formation of galaxies in a ACDM uni- 
verse being conducted by the Virgo consortiimi - the Galaxies- 
Intergalactic Medium Interaction Calculation or GiMiC. The goal 
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Figure 12. Median baryon fraction, /t, = {My^ / M200) / (f^b/f^m), of the simulated lialoes. Tlie dotted line indicates the 90 per cent level expected in the 
non-radiative regime; the arrows indicate the wind launch velocity. Box and whisker elements show the median and 10**", 25**", 75*'' and 90"^ percentiles of 
the total. Left: resolution test, comparing the two realisations of the —2a region. Agreement between the high- and intermediate-resolution runs for the total 
baryon fraction is achieved at a circular velocity that decreases with redshift; at z = this scale is ~ 50 km s^^. Right: coloured curves decompose the 
baryon fraction into the cumulative contribution made by stars (dark blue), equation of state gas {green), cold gas (cyan), warm gas (orange) and hot gas (red), 
in that order. Below 160 km s"'^, results are drawn from the high-resolution —2(7 region; above this scale they are drawn from all five intermediate-resolution 
regions. Dotted sections of the curves (Vc < 50 km s~^) indicate the regime where we believe the high-resolution simulations to be unreliable. 



of the programme is to investigate the joint evolution of galax- 
ies and the intergalactic medium (IGM). The simulations follow 
the evolution of nearly spherical regions of radius ~ 20 h^^ Mpc, 
drawn from the Millennium Simulation (Springel et al. 2005). The 
regions have mean overdensities of (—2, —1, 0, +1, +2)a, where 
a is the root mean square overdensity on the scale of the regions 
at z = 1.5. The rest of the Millennium Simulation volume is sim- 
ulated at much lower resolution using coUisionless particles. This 
provides the correct tidal forces on the high-resolution regions. 

The five GiMIC regions were followed to 2; = at inter- 
mediate resolution (rrigas ~ 1.16 x 10^ M0). Of the high- 
resolution simulations (mgas = 1.45 x 10^ Mq), the —2a 
region was continued to z = 0, whilst the (— 1,0,+1)(t regions 
were stopped at 2 — 2. In the highest resolution simulations the 
Jeans mass in the IGM after reionisation is marginally resolved and 
thus these simulations can, in principle, track the formation of all 
galaxies whose cooling is dominated by Hi. 

Our simulation strategy has a number of advantages: (i) it sam- 
ples the entire range of large-scale environments - including rare 
objects such as massive cluster haloes and voids - allowing us to ex- 
amine the dependence of galaxy properties on large-scale overden- 
sity; (ii) it allows accurate integration to 2; = since fluctuations 
on the scale of the computational volume remain linear; (iii) by 
suitably averaging over the five regions, it allows an estimate of the 
statistical properties of representative cosmological volumes. 

The simulations were carred out using the GADGETS code, a 
substantial upgrade of GADGET2 (Springel 2005) that includes; 

(i) a recipe for star formation designed to enforce a local 
Kennicutt-Schmidt law (Schaye & Dalla Vecchia 2008); 



(ii) stellar evolution and the associated delayed release of 11 
chemical elements (Wiersma et al. 2009b); 

(iii) the contribution of metals to the cooling of gas, computed 
element-by-element, in the presence of an imposed UV-background 
(Wiersma et al. 2009a); 

(iv) galactic winds that pollute the IGM with metals and can 
quench star formation in low-mass haloes (Dalla Vecchia & Schaye 
2008). 

For these simulations we do not, however, follow the evolution of 
black holes or feedback effects associated with them. 

The simulations produce the correct number density of galax- 
ies with mass greater than about 10® Mq . At 2; = 2 the stellar mass 
function is consistent with observations. At 2; = there are too 
many galaxies of very low and very high mass and not enough of 
intermediate mass. The excess at the massive end reflects the ab- 
sence of a heating mechanism, such as that often ascribed to AGN 
feedback, to prevent too much gas cooling. At small and interme- 
diate masses the discrepancy with observations reflects inadequa- 
cies in our treatment of galactic winds. For example, the paucity of 
galaxies with stellar mass ~ lO'^" Mq seems to be directly related 
to our simple wind model. 

In this paper we have concentrated on the star formation prop- 
erties of the simulations. These are governed by an interplay be- 
tween accretion, gas cooling, photoheating by the imposed ion- 
ising background and galactic winds. The specific star formation 
rate in our simulations is converged for z < 6 and we estimate 
that the value of p*{z) in the high-resolution simulations is re- 
liable for z < Zrcion = 9. Indeed, the reionisation process is 
"visible" in these simulations as a depression of the star formation 
rate in low mass haloes. In all five GiMIC regions, Pi,{z) increases 
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slowly with time, reaches a broad maximum around 2; ~ 2, and 
then falls rapidly by a factor of ~ 6 to the present day. In general, 
the global value of p*(«), obtained by a weighted average over the 
individual regions, is in reasonable agreement with the observa- 
tional compilation of Hopkins (2007). However, the decline in 
at low-redshift in the simulation is not as pronounced as in the ob- 
servations, largely because of excessive star formation in massive 
haloes. 

One of the main results of this paper is the strong environmen- 
tal dependence we find of the star formation rate density, p*(2), 
which is manifest as a large variation amongst the five GiMic re- 
gions. Although the shapes of the p* {z) curves are similar for all re- 
gions, the offset in amplitude between the —2a to the +2a region is 
approximately an order of magnitude. The Oct region closely tracks 
the global value reconstructed from a weighted average of all five 
regions and is typically a factor of 2 below the +2a region. Even 
when normalised by total mass rather than volume, p* still varies 
by more than a factor of 3 between the two extreme regions. This 
environmental dependence makes it difficult to obtain reliable ob- 
servational estimates, particularly at high-redshift. For example, at 
a = 6 the 18 Mpc GiMic spheres correspond to 320 arcmin^ 
on the sky, which is a considerably larger area than the (5 — 20 
arcmin^) covered by the Hubble Ultra Deep Field. 

At a given redshift, we find that the specific star formation 
rate in dark matter haloes, M,/M2oo (the "halo sSFR"), depends 
only on halo mass, not on large-scale environment. At all redshifts, 
the specific star formation is most efficient in haloes of circular 
velocity u,; = 200 - 250 km Above this value, the sSFR 
varies only little with Vc, but below this value it drops rapidly due 
to the action of galactic winds that deplete these haloes of baryons 
to below ~ 10 per cent of the cosmic mean by z = 0. Above this 
critical value of Vc, the sSFR is well described by the spherically 
symmetric cooling model of White & Frenk (1991). Because of 
their generally higher baryon fractions and gas densities, haloes of 
a given mass tend to have higher sSFR at earlier times. 

Since the halo sSFR is independent of environment, the rea- 
son for the strong dependence of the star formation rate density 
with environment is a strong dependence of the halo mass func- 
tion on enviroimient. There are always more haloes and galax- 
ies in the +2(t region than in the —2<r region. For example, the 
volume-normalised number density of galaxies with stellar mass 
> 10^° ft"^ M© at 2; = is a factor of ~ 3 greater in the 
+2cr than the — 2ct region. In addition, the most massive haloes 
are strongly biased towards regions of high overdensity. Since the 
sSFR in our simulations is always dominated by high mass haloes 
with Vc = 200 — 250 km s^^, even when normalised to the total 
mass in each region, the +2a region always has a higher specific 
star formation rate density than the —2(t region. The combined ef- 
fect of halo number density and star formation efficiency results 
in star formation at the present time that is dominated by galac- 
tic haloes (A/2oo ~ 10^^ h^^ M©) in the void region, by groups 
(M200 ~ 10^^"^'' M0) in the intermediate density regions, 
and by clusters, M200 ~ 10^® M©, in the highest density re- 
gion. 

The plateau in the global p* (z) relation at 6 > z > 2 is due to 
the balance between the advancing formation of dark matter haloes 
and the decline in the sSFR in haloes of a given mass (see Figs. 8 
and 11). Eventually, the latter dominates and p* plummets in all 
regions. The lack of feedback effects associated with the growth 
of central black holes may affect the quantitative, but perhaps not 
the quaUtative, behaviour of these trends. This view is supported 
by the agreement, to within a factor of ~ 2, between the overall 



star formation rate densities inferred for the Millennium Simulation 
from our simulations and from the semi-analytic models of Bower 
et al. (2006); Croton et al. (2006) and De Lucia et al. (2006), which 
do include the effects of AGN feedback. 

In common with the conclusions of recent studies employing 
semi-analytic modelling, we find that the growth of galaxies in our 
simulations is segregated by stellar mass: low-mass galaxies con- 
tinue to grow at all epochs due to ongoing star-formation, whilst 
more massive galaxies have largely concluded their star formation 
at 2: > 1. This apparent "downsizing" of the star formation activity 
occurs in spite of the absence of AGN feedback in the simulations. 

The star formation rate density in the simulations is always 
dominated by the progenitors of galaxies that today reside in 
the most massive haloes (group/cluster haloes with M ^ 5 x 
10^^ M.q). Yet, these stars were predominantly formed either 
in dwarf haloes (M < 5 x lO" h'^ Mq) before z ~ 3.5 or 
in galaxy-sized haloes after that. Galactic winds are an effective 
mechanism for quenching star formation to the observed levels in 
these smaller systems. The integrated stellar mass density at 2 = 
agrees with observational estimates to within a factor of < 2. 

Some properties of the galaxies in the simulations can be com- 
pared with observations. Galaxies form first in the high density re- 
gion. Between z ~ 9 — 6, they are very compact (radii ^ ih^^ 
kpc), massive (M* ~ 10^^ h^^ Mq), contain relatively old stars 
(~ 0.2 Gyr old), have high metallicity {Z ~ Z©) and are over- 
abundant in a-elements. These galaxies are reminiscent of those 
recently found in the GOODS fields by Wiklind et al. (2008). 

In summary, the GiMIC simulations provide a useful resource 
to investigate the evolution of and interaction between galaxies and 
the intergalactic gas in different large-scale environments. In this 
first paper of a series, we have explored the star formation proper- 
ties of haloes and some of the properties of the galaxies that form 
within them. In subsequent papers, we will examine further prop- 
erties of both galaxies and intergalactic gas. 
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APPENDIX A: METHODS IN DETAIL 

In this appendix we expand on the technical details of the genera- 
tion of the initial conditions for GiMIC simulations. Secondly, we 
discuss how we address the existence of an artificial boundary be- 
tween the high-resolution, baryon-fiUed region of each simulation, 
and the extemal low-resolution region that is simulated with colUs- 
sionless dynamics. As part of this latter discussion, we describe our 
technique for combining results from the five GiMIC simulations to 
construct estimates of properties of the Millennium Simulation as 
a whole. 

Al Generation of the initial conditions 

Here we provide a more detailed description of the generation of 
the initial conditions, than the brief overview that was given in Sec- 
tion 2.2. As described there, we aim to trace a representative sample 
of the (500 h^^ Mpc)"^ Millennium Simulation volume with five 
separate simulations, and also to follow the evolution of rare struc- 
tures such as voids and clusters that are only found in such large pe- 
riodic volumes. Clearly, it is desirable to perform these simulations 
at high resolution, and we aimed to have a sufficiently low parti- 
cle mass that the Jeans scale in the intergalactic medium (after the 
epoch of reionisation) would be resolved in our "high-resolution" 
runs (which have 8 times more particles than our 'intermediate- 
resolution' rans). In practice, this limits the size of the regions to 
~ 18 hr^ Mpc. However, our additional constraint that the +2(7 
region should be centred on a rich cluster at 2: = required the size 
of this region to be increased to 25 h^^ Mpc; this is because the 
mass of a sphere of mean density and radius 18 Mpc at 2; = 
is 1.7 X 10^® Mq, a mass comparable to that of the largest 
cluster in the Millennium Simulation. We describe the generation 
of the initial conditions in two parts: i) the selection of the regions 
to be resimulated, and ii) the generation of the particle distribution 
and the displacement field. 
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Al.l Selecting the five regions 

The four (—2, —1, 0, +1)(t spheres were selected from a target list 
of 10^ randomly placed spheres of radius 18 h^^ Mpc in the Mil- 
lennium Simulation volume at z = 1.5. Since hy z = 1.5 the 
overdensity distribution has become significantly skewed (i.e. non- 
Gaussian), we defined the overdensities of the (—2, —1, 0, +l)a 
regions as those that correspond in terms of the rank ordering by 
overdensity, to a Gaussian distribution; the overdensities for the 
(-2,-1,0, +l)cr regions are (-0.407,-0.243,-0.032,0.241) 
respectively. From the complete list of 10 ' spheres a candidate list 
of centres was generated for each of the regions by selecting all 
spheres with overdensities within 0.002 of the respective Umits. 

The procedure for the +2tr sphere was different because of 
the requirement that a rich cluster should form at its centre by 
z = 0. The starting point was to identify the particles that form the 
126 most massive groups (that have a characteristic separation of 
~ 100 Mpc) in the Millennium Simulation volume at ^ = 0, 
identified by a friends-of-friends group finder (Davis et al. 1985) 
with linking length b = 0.2. These particles were traced back to 
z = 1.5 and their centres of mass used as centres for spheres of 
radius 25 Mpc. The +2(7 sphere was then selected from a ran- 
domly sorted list of twelve candidates that had an overdensity close 
to the appropriate value. Since the first cluster had a mass at z = 
of Af2()() = 4 X lO^'' h^^ Mq, which is rather low for a "rich" 
cluster, we chose the second sphere on the list, with a mass of 
M200 = 8 X 10" M© at = 0. This cluster is relatively 
isolated, with no material in the 25 Mpc belonging to another 
rich cluster. 

Having selected the centre for the +2(t sphere, and having 
made candidate lists of centres for the four other spheres, we se- 
lected our final five sphere centres from a list of all five centres 
generated as follows. We selected a sphere, at random, from each of 
the candidate lists of centres for the (— 2, — 1, 0, +1)(7 regions and 
added the chosen +2a sphere centre to form a quintet of candidate 
centres. Each quintet was considered eligible for the final list if it 
satisfied two criteria: (i) none of the five centres had one or more 
of their x,y or z coordinates within 50 h^^ Mpc of the periodic 
boundary of the simulation volume and (ii) each centre was at least 
200 Mpc away from the other four sphere centres in the quin- 
tet. The first criterion was imposed to maintain the same coordinate 
system as the Millennium Simulation, and also allow the use of an 
intermediate-resolution particle-mesh (PM) force algorithm imple- 
mented within our simulation code, which requires that the isolated 
mesh used for this calculation should not cross the periodic bound- 
aries of the simulation volume. The second criterion was chosen to 
prevent two or more spheres from being near neighbours, and to 
force the five spheres to extend over a significant range of the Mil- 
lennium Simulation voliune, so as to be relatively independent of 
each other. It is possible to find many quintets satisfying these con- 
ditions; we chose the first quintet generated by the code that used 
a pseudo-random number generator to search for candidates. The 
centres and radii of all five spheres are presented in Table 1. 

In order to maximise the load and memory efficiency of the 
intermediate-range PM algorithm in our code, which automati- 
cally resizes to encompass all high-resolution particles, we chose to 
make five separate sets of initial conditions, one for each of the re- 
gions chosen; this is because the inclusion of all five spheres within 
a single set of initial conditions would create an isolated mesh fill- 
ing most of the simulation volume. Such a mesh would incur a sig- 
nificant memory cost and offer little speed up relative to the main 
PM algorithm. 



A1.2 Generating the particle load and applying the 
displacement field 

The procedure for making the initial conditions is identical for all 
five regions and consists of two main stages: (i) the creation of a 
multi-mass particle distribution to represent the uniform, unper- 
turbed particle distribution, and (ii) the recreation of the displace- 
ment field used by the Millennium Simulation, with the addition of 
extra short wavelength power in the resimulated regions, sampled 
from the same power spectrum, down to the appropriate Nyquist 
frequency. The multimass uniform particle distribution is based 
upon a cubic mesh of 320'' cells that encloses the entire simulation 
volume. Each cell of the mesh is treated independently. Particles, 
which at this point are composite and represent both dark matter 
and gas, are placed in each cell so that the mean density of all cells 
is equal to the matter density of the Universe. A total of n'^ equal 
mass particles is placed in each cell, arranged as a cubic grid with a 
spacing of 1/n of the cell size and with the centre of mass of all the 
particles coinciding with the cell centre. The value of n varied from 
32, for the high-resolution region encompassing the sphere itself, to 
unity for cells that are distant from the sphere and are needed only 
to provide the correct tidal forces on the region of interest. 

The value of n, as a function of cell position, was determined 
by the following algorithm: all of the particles from the Millennium 
Simulation inside or within 0.5 h^^ Mpc (1 h^^ Mpc for the +2cr 
sphere) of a particular GiMic sphere atz = 1.5 were traced back to 
the grid cells from which they originated at very high redshift. Any 
cell occupied by one or more of these particles is labelled as a high- 
resolution cell. A check is made to determine if the high-resolution 
cells form a single contiguous region, where contiguous cells are 
defined as those sharing a face. If there were non-contiguous cells 
present, then any cells contiguous to a high-resolution cell are addi- 
tionally defined as high-resolution cells, and this process is contin- 
ued until all high-resolution cells are contiguous. A set of boundary 
layers is then defined as follows: level 1 boundary cells are defined 
as those cells adjacent (i.e. sharing a face, an edge or a vertex) to 
a high-resolution cell. Level 2 boundary cells are then adjacent to 
level 1 boundary cells, and so on until level 9, which in addition 
includes all cells that have not yet been assigned. 

We then realise the zoomed initial conditions at two reso- 
lutions; since the original Millennium Simulation represents our 
low-resolution basis, we refer to the two GiMIC realisations as 
high- and intermediate-resolution. For the high-resolution reali- 
sation, we used n = 32 in the high-resolution cells, while the 
neighbouring boundary cells had n = (20, 15, 10, 8, 6, 4, 3, 2, 1) 
for the boundary levels 1-9 respectively. This yields a mass for 
the composite (i.e. dark matter and gas) particles in the high- 
resolution region of 8.08 x 10^ Mq. For the intermediate- 
resolution realisation, we effectively halve the value of n in each 
level (odd values are rounded up), yielding initial conditions with 
composite particle masses a factor of eight greater than the high- 
resolution case, 6.46 x lO'^ Mq. We therefore resolve haloes 
of mass 10^^ Mq (comparable to the Milky Way's halo) with 
~ 19, 000 particles at intermediate-resolution, and ~ 150, 000 at 
high-resolution. The total number of particles used in each simula- 
tion is presented in Table 1 . 

Having created the particle distribution, a displacement field 
was calculated for each particle using the techniques for making 
resimulation initial conditions described in Power et al. (2003) and 
Navarro et al. (2004). The displacement consists of two parts; on 
large scales the same mode amplitudes and phases were used as in 
the Millennium Simulation, whilst short wavelength power down 
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to the Nyquist frequency was added in a cubic region of dimension 
50 — 76 Mpc containing all of the high-resolution cells. The 
break between the long wavelength and short wavelength power 
was chosen such that only wavenumbers k — {kx,ky,k:,) satis- 
fying max(|A:2,|, \ky\, |fcz|) < nh Mpc~^ were retained from the 
Millennium Simulation. This scale is sufficient to ensure that haloes 
resolved with at least a few hundred particles in the Millennium 
Simulation are present in the GiMIC resimulation. A large mesh of 
2048^ cells was used for all Fourier transforms to minimise inter- 
polation errors. 

Finally, the high-resolution composite particles were split into 
gas and dark matter particles, with masses appropriate to yield a 
cosmic baryon fraction of fib/i^m = 0.045/0.25 = 18 per cent. 
The gas and dark matter particles were offset in each dimension 
from their common centre of mass by a mass-weighted fraction of 
one interparticle separation. 

A2 Sampling & combining regions 

The use of zoomed initial conditions complicates the interpreta- 
tion of simulations because of the presence of an artificial boundary 
around the high-resolution region. Inside it, baryons are modelled 
explicitly with SPH particles; outside it, dark matter and baryons 
are represented by composite particles whose dynamics are colli- 
sionless. Clearly, the coUisionless region should be excluded from 
any analysis, but the existence of the boundary has other indirect 
effects. For example, the absence of gas beyond the boundary arti- 
ficially reduces the mass of gas that might accrete onto peripheral 
haloes. In addition, the lack of external pressure allows the easy es- 
cape of galactic outflows into the vacuum of the boundary region. 

In the analyses presented in this paper that measure proper- 
ties on a particle-by-particle basis, we consider only those particles 
that are within 18 Mpc (comoving; 25 Mpc for the +2(j 
region) of the centre of mass of all baryonic and high-resolution 
dark matter particles; in the case of analyses where we consider the 
properties of haloes or galaxies, we consider those whose centre-of- 
mass is within this volume. For 2 > 1.5 this automatically leaves a 
small boundary layer of SPH and high-resolution dark matter par- 
ticles external to the sphere for "safety", since a padding region 
was included in the initial conditions by design. However, at lower 
redshifts the morphology and comoving volume of each region de- 
viates slightly from that of the sphere that defined them at z = 1.5. 
This is due to the varying overdensity of each region, since the 
underdense regions expand and the overdense regions contract in 
comoving space. However, we have explicitly checked that haloes 
towards the edge of the volume do not differ significantly from 
those at the centre, for instance in terms of their baryon fraction, 
star formation rate or metallicity, and so we do not consider this an 
issue for the analyses we present. In addition, the use of a simple 
and time-invariant sampling geometry (a sphere of constant comov- 
ing radius) allows various quantities, such as the star formation rate 
density p*, to be normalised by volume in an unambiguous fashion. 

In Fig. Al we present the overdensity evolution of the spheri- 
cal regions within each simulation. At early times, the regions nec- 
essarily have the mean density of the Universe, but their enclosed 
overdensities clearly diverge as they evolve. This figure highlights 
the advantage of zoomed initial conditions as a means to probe 
varying environments, since a periodic volume maintains 5 = 
at all epochs. The dip in the overdensity of the +1(7 region at 
a ~ 0.5 {z ~ 1) results from the departure of a high-mass halo 
with a large peculiar velocity from the spherical "analysis" region. 

At 2 = 1.5 the distribution of overdensities over spherical 
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Figure Al. The overdensity evolution of the five spherical regions at in- 
termediate resolution. The vertical dotted line denotes the epoch at which 
the regions were selected (z = 1.5). Notice the drop, after z = 1.5, in 
overdensity of the +lcr region; this results from a massive halo with high 
peculiar velocity leaving the spherical region. 
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Table Al. Numerical weights applied to each GiMIC region when extrapo- 
lating statistics to the whole (500 Mpc)'' Millennium Simulation vol- 
ume. The net factor (final column) is the product of the summation weight 
assigned to each region (first column) and the volume weighting (second 
column). The latter is applied to account for the larger volume of the +2(7 
region. 



regions of size comparable to the GIMIC regions is close to Gaus- 
sian. Ideally, to obtain the universal average of a quantity which is a 
function of overdensity, it is necessary to integrate over overdensity 
with a Gaussian weighting. Because the GiMIC simulation provides 
only 5 overdensities. (— 2, — 1, 0, +1, +2)(7. we approximate the 
integral as a sum and choose the weights for the five regions so that 
the summation correctly integrates any odd polynomial of overden- 
sity and even polynomials up to and including a quartic power in 
overdensity. These weights are given in the second column of Ta- 
ble Al. An additional complication is that the +2(j region has a 
different size, so an additional weight is needed to compensate for 
the larger volume (Table Al, third column). The actual weights are 
given in the fourth column Table Al . 

These weights can only be approximate because the over- 
density distributions are not exactly Gaussian, and strictly the re- 
gions only correspond to (— 2, — 1, 0, +1, +2)(j deviations from 
the mean a\.z = 1.5. In Fig. A2, we show the halo multiplicity func- 
tion, defined as the fraction of all mass bound into haloes of mass in 
the interval (log M, log M + d log M). This description of the halo 
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Figure A2. The halo multiplicity function in the five regions at intermedi- 
ate resolution. The black histogram is the halo multiplicity function yielded 
by our weighting procedure used to generate estimates for the Millennium 
Simulation volume from the GiMIC regions. For comparison, we show the 
multiplicity function derived from the Reed et al. (2007) mass function (dot- 
ted line). 




population is attractive since it has a smaller dynamic range than the 
halo mass function, and so more clearly highlights the differences 
between the halo populations of the GiMIC regions. Inspection of 
this plot shows that our weighted average (black histogram) con- 
structed from the individual regions (coloured histograms) recovers 
the halo multiplicity function (and therefore the halo mass function 
also) of the Millennium Simulation (dotted line) reasonably well 
at all times. It is this fact that ultimately justifies our weighting 
scheme. 



